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ABSTRACT 

In this thesis we present a numerical study of general relativistic boson stars in both spherical 
symmetry and axisymmetry. We consider both time-independent problems, involving the solution 
of equilibrium equations for rotating boson stars, and time-dependent problems, focusing on black 
hole critical behaviour associated with boson stars. 

Boson stars are localized solutions of the equations governing a massive complex scalar field 
coupled to the gravitational field. They can be simulated using more straightforward numerical 
techniques than are required for fluid stars. In particular, the evolution of smooth initial data for 
a scalar field tends to stay smooth, in sharp contrast to hydrodynamical evolution, which tends to 
develop discontinuities, even from smooth initial conditions. At the same time, relativistic boson 
stars share many of the same features with respect to the strong-field gravitational interaction as 
their fermionic counterparts. A detailed study of their dynamics can thus potentially lead to a 
better understanding of the dynamics of compact fermionic stars (such as neutron stars), while the 
relative ease with which they can be treated numerically makes them ideal for use in theoretical 
studies of strong gravity. 

In this last vein, the study of the critical phenomena that arise at the threshold of black hole 
formation has been a subject of intense interest among relativists and applied mathematicians over 
the past decade. Type I critical phenomena, in which the black hole mass jumps discontinuously 
at threshold, were previously observed in the dynamics of spherically symmetric boson stars by 
Hawley and Choptuik We extend this work and show that, contrary to previous claims, the 

subcritical end-state is well described by a stable boson star executing a large amplitude oscillation 
with a frequency in good agreement with that predicted for the fundamental normal mode of the 
end-state star from linear perturbation theory. 

We then extend our studies of critical phenomena to the axisymmetric case, studying two dis- 
tinct classes of parametrized families of initial data whose evolution generates families of spacetimes 
that "interpolate" between those than contain a black hole and those that do not. In both cases 
we find strong evidence for a Type I transition at threshold, and are able to demonstrate scaling of 
the lifetime for near-critical configurations of the type expected for such a transition. This is the 
first time that Type I critical solutions have been simulated in axisymmetry (all previous general 
relativistic calculations of this sort imposed spherical symmetry). 

In addition, we develop an efficient algorithm for constructing equilibrium configurations of 
rotating boson stars, which are characterized by discrete values of an angular momentum parameter, 
k (an azimuthal quantum number). We construct families of solutions for k — 1 and k — 2, and 
demonstrate the existence of a maximum mass in each case. 
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CHAPTER 1 
INTRODUCTION 

This thesis is concerned with the numerical simulation of boson stars within the framework of 
Einstein's theory of general relativity. Boson stars are self-gravitating compact objects 1 composed 
of scalar particles |3j , and are the bosonic counterparts of the more well-known compact fcrmionic 
stars, which include white dwarfs and especially neutron stars. In contrast to fermionic stars, 
there is no observational evidence that boson stars exist in nature, nor has the existence of any 
fundamental scalar particle yet been verified by experiments. It has been proposed, however, that 
boson stars could be a candidate for, or at least make up a considerable fraction of, the dark matter 
in our universe 0]. Were this true, it is quite probable that studies of their properties would lead 
to a better understanding of astrophysical phenomena. However, at the current time boson stars 
remain purely theoretical entities. 

Their current hypothetical nature notwithstanding, boson stars might indeed exist in our uni- 
verse, and, more importantly for this thesis, they are excellent matter models for the numerical 
study of compact objects in strong gravitational fields, a subject that continues to be a key con- 
cern of numerical relativity 5 4 (numerical relativity: the numerical simulation of the Einstein field 
equations, as well as the field equations for any matter fields being modeled). The boson stars 
we consider can be viewed as zero-temperature, ground-state, Bose-Einstein condensates, with 
enormous occupation numbers, so that the stellar material is described by a single complex scalar 
field, <f>(t,x), that satisfies a simple, classical Klein-Gordon equation. Many of the nice modeling 
properties associated with boson stars derives from the fact that the dynamics is governed by 
a partial differential equation (PDE) that does not tend to develop discontinuities from smooth 
initial data. This is not the case for fermionic stars, which are usually treated as perfect fluids, 
and which thus satisfy hydrodynamical equations with phenomenologically-determined equations 
of state. Relativistic fluid evolution will generically produce shock waves and other discontinuities, 
even from smooth initial data. Therefore, today's state-of-the-art relativistic hydrodynamics codes 
use sophisticated high resolution shock capturing (HRSC) schemes in order to more accurately and 
efficiently simulate the fluid physics. The fact that the solutions of the Klein-Gordon equation stay 
smooth is a tremendous boon when it comes to discretizing the equation in a stable and accurate 
manner. 

In addition, the effective integration of HRSC schemes with other advanced numerical techniques 
such as adaptive mesh refinement [§] (AMR) is still in its infancy, at least for the case of general 
relativistic applications. For bosonic matter, on the other hand, it has proven relatively easy 
to use rather generic AMR algorithms (such as that of Berger and Oliger |H]) in conjunction 
with straightforward second-order finite-difference discretization to achieve essentially unbounded 

1 By compact, we mean gravitationally compact, so that the size, R, of the star is comparable to its Schwarzschild 
radius, Rg. The Schwarzschild radius associated with a mass M is Rg = (2G/c 2 )M, where G is Newton's gravita- 
tional constant, and c is the speed of light. 
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dynamical range in an efficient fashion 0EHH]. 

Of course, we do not primarily study boson stars because they are easy to simulate. As men- 
tioned already, a good fraction of the work in numerical relativity is concerned with the strong- field 
dynamics of gravitationally compact objects @JE3E]- Despite the large amount of effort that 
has been devoted to this subject, it is fair to say that much remains to be learned, and much of 
what remains to be discovered is likely to be found through simulation. Although the strong-field 
gravitational physics of boson stars may not compare in detail with that of fermionic stars in all 
respects, there are clearly some key features of the usual stars that are shared by their bosonic 
counterparts. For example, in analogy with relativistic fermionic stars, spherically symmetric bo- 
son stars typically come in one-parameter families, where the parameter can be viewed as the 
central density of the star (or an analogue thereof). Moreover, in relativistic cases the boson star 
families share with the fermionic sequences the surprising property that there comes a point when 
increasing the density of the star at the center actually decreases the total gravitating mass of the 
star. In both cases this leads to maximal masses for any given family of stars. In addition, in both 
instances, stars near that limit are naturally very strongly self-gravitating. 

Suffice it to say, then, that a careful study of general relativistic boson stars is likely to lead 
to insights into strong-field gravitational physics, even if there are no immediate astrophysical 
applications, and that this is the primary motivation for the calculations described below. In 
particular, and again as with the fermionic case, black hole formation is a crucial process that can, 
and will, generically occur in the strong-field dynamics of one or more boson stars. For the most 
part, this is simply because (1) as gravitationally compact objects, boson stars are, by definition, 
close to the point of collapse, and (2) the process of gravitational collapse involving matter with 
positive energy tends to be very unstable in many senses, including the fact that black hole areas 
can only increase. 

Over the past decade or so, the careful study of gravitational collapse and black hole formation 
has lead to the discovery of black hole critical phenomena 0, wherein the process of black hole 
formation, studied in solution space, takes on many of the features of a phase transition in a 
statistical mechanical system. A primary goal of the work presented in this thesis is to study so- 
called Type I critical phenomena of boson stars in both spherical symmetry and axisymmetry, and 
this will be explained in more detail below. The simulations are carried out via finite difference 
solution of the governing PDEs. We also develop and apply a new algorithm for constructing 
solutions of the time-independent form of the PDEs in axisymmetry; these solutions represent 
relativistic rotating boson stars. 

1.1 An Overview of Boson Stars 

The study of boson stars can be traced back to the work of Wheeler. Wheeler studied self- 
gravitating objects whose constituent element is the electromagnetic field and named the resulting 
"photonic" configurations geons |12| . Wheeler's original intent was to construct a self-consistent, 
classical and field-theoretical notion of body, thus providing a divergence-free model for the New- 
tonian concept of body. In the late 1960s, Kaup 3 adopted the geon idea, but coupled a massive 
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complex scalar field, rather than the electromagnetic field, to general-relativistic gravity. Assuming 
time-independence and spherical symmetry he found solutions of the coupled equations which he 
called Klein-Gordon geons. Subsequently, Ruffini & Bonazzola |KS| studied field quantization of a 
real scalar field and considered the ground state configurations of a system of such particles. The 
expectation value of the field operators gives the same energy-momentum tensor as those given 
by Kaup, and hence the different approaches followed in the two studies give essentially the same 
macroscopic results. 

Later, these Klein-Gordon geons were given the name boson stars, and the nomenclature boson 
star now generally refers to compact self-gravitating objects that are regular everywhere and that 
are made up of scalar fields. Variations of the original model studied by Kaup and Ruffini & 
Bonazzola include self-interacting boson stars (described by a Klein-Gordon field with one or more 
self-interaction terms) , charged boson stars (boson stars coupled to the electromagnetic field) and 
rotating boson stars (boson stars possessing angular momentum), to name a few. When stable, all 
of these objects are held together by the balance between the attractive gravitational force and a 
pressure that can be viewed as arising from Heisenberg's uncertainty principle, as well as any explicit 
repulsive self-interaction between the bosons that is incorporated in the model. Depending on the 
mass of the constituent particles, and on the value of the self-interaction coupling constant(s), the 
size of the stars can in principle vary from the atomic scale to an astrophysically- relevant scale ^Oj • 
As mentioned previously, any given boson star is typically only one member of a continuous family 
of equilibrium solutions. In spherical symmetry, the family can be conveniently parametrized using 
the central value of the modulus of the scalar field, in analogy to the central pressure of perfect 
fluid stars. Moreover, the equilibrium configurations generally have an exponentially decaying tail 
at large distances from the stellar core, in contrast to fluid-models stars which tend to have sharp, 
well-defined edges. 

1.1.1 Stability of Boson Stars 

The dynamical stability of equilibrium, compact fermionic (fluid) stars against gravitational col- 
lapse can be studied using the linear perturbation analysis of infinitesimal radial oscillations that 
conserve the total particle number, N, and mass/energy, M. One important theorem in this re- 
gard concerns the transition between stable and unstable equilibrium PP-305]. The theorem 
states that a perfect fluid star with constant chemical composition and constant entropy per nucleon 
becomes unstable with respect to some radial mode only at central densities p(0) such that 

dM {pi 0),s,..) = Q) 



d P (p) 

dN(p(0),s,---) 

d P (o) 



0. (1.2) 



In other words, a change in stability can only occur at those points in a curve of total mass M(p(0)) 
vs p(0), where an extrema is attained. This results directly from the fact that the eigenvalues, c 2 , 
of the associated pulsation equation change sign at those points. 

Similar results hold for various boson star models. As just mentioned, in the case of boson stars 



CHAPTER 1. INTRODUCTION 



4 



we use the central value of the scalar field modulus, denoted 4>o(0), as the parameter for the family 
of solutions, and it has been shown that the pulsation equation has a zero mode at the stationary 
points in the M(<fio(0)) plot ^JEl- Numerical verification of the instability of configurations 
past the mass maximum using the full dynamical equations in spherical symmetry was studied in 
PUDS] and E]. 

1.1.2 Maximum Mass of Boson Stars 

As already mentioned, another important feature of stellar structure that is largely due to relativis- 
tic effects is the existence of a maximum allowable mass for a particular species of stars. Depending 
on the mechanism of stellar pressure (degenerate electron pressure for white dwarfs, degenerate neu- 
tron pressure for neutron stars, Heisenberg uncertainty principle for boson stars) these values are 
different. White dwarfs and neutron stars share the same dependence on the mass of the con- 
stituent particles (~ M^/m 2 ), while the dependence for boson stars is quite different (~ M 2 j/to). 
The origin of the difference can be understood via the following heuristic argument. The ground 
state stationary boson stars are macroscopic quantum states of cold, degenerate bosons, whose 
existence is the result of the balance between the attractive gravitational force and the dispersive 
nature of the wave function. By the uncertainty principle, if the bosons are confined in a region of 
size i?, we have pR ~ h, where p is the typical momentum of the bosons. For a moderately rela- 
tivistic boson we have p ~ mc and hence R ~ H/ {mc). Equating this with the Schwarzschild radius 
we have H/(mc) ~ 2GM max /c 2 . Hence M max ~ (Hc)/(2Gm) = O.bM^/m. Numerical calculation 
shows that M max « 0.63iM 2 Jm (see Sec. 14.3. 4"t . surprisingly close to the above estimate. 

The maximum mass of neutron stars can be estimated in a similar way. 2 The existence of 
these stars is the result of the balance between the attractive gravitational force and the pressure 
due to degenerate neutrons (fermions). Suppose there are N fermions confined in a region of size 
R. Then by Pauli's exclusion principle, each particle occupies a volume 1/n, where n = N/R 3 
is the number density. Effectively, each particle has a size of R/N 1 / 3 . Again, by the uncertainty 
principle we have pR/N 1 ^ 3 ~ h. Following the same argument as for the boson star case, we have 
R ~ hN^ 3 /(mc), and hence 2GM max /c 2 - HN 1 ' 3 /^) ~ Mf max 1/3 /(™ 4/3 c). Thus we have 
A^max ~ 0.35Mpi 3 /m 2 . In contrast to the bosonic case, then, the maximum mass of fermionic stars 
scales as M max ~ M v 3 /m 2 . 

1.1.3 Rotating Boson Stars 

Although spherically symmetric boson star solutions were found as early as 1968, for many years 
it was unclear whether solutions describing time-independent boson stars with angular momentum 
existed or not. The first attempt to construct such stars was due to Kobayashi, Kasai & Futamase 
in 1994 These authors followed the same approach as Hartle and others E^IEHI, m which the 
slow rotation of a general relativistic star is treated via perturbation of a spherically symmetric 
equilibrium configuration. In their study they found that slowly rotating boson stars solutions 

2 1201 gives another heuristic argument due to Landau (1932), which applies to both neutron stars and white 
dwarfs. The current argument applies to neutron stars only; for white dwarfs the mass of the star is dominated by 
baryons, while the pressure is provided by electrons. 
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coupled with a U(l) gauge field (charged boson stars) do not exist, at least perturbatively. Shortly 
thereafter, however, it was demonstrated that rotating boson stars could be constructed on the 
basis of an ansatz which leads to a quantized (or perhaps more properly, discretized) angular 
momentum. It was later understood that rotating boson stars have quantized angular momenta 
and that the concomitantly discrete nature of the on-axis regularity conditions prohibit a continuous 
(perturbative) change from non-rotating to rotating configurations. 

A year later, Silveira & de Sousa |23], following the approach of Ferrell & Gleiser j^S], succeeded 
in obtaining equilibrium solutions of rotating boson stars within the framework of Newtonian 
gravity. Specifically, they adopted the ansatz 3 

<t>(t,r,9,cp) =Mr,6)e i(ut+k ' p) , (1.3) 

where k is an integer (we use the symbol k, instead of the symbol to, that is commonly used in 
quantum mechanics, to avoid confusion with the particle mass, to), so that the stress-energy tensor 
Tt v is independent of both time t and the azimuthal angle ip. In the non-relativistic limit 4 the 
governing field equations constitute a coupled Poisson-Schrodinger (PS) system (for simplicity, we 
have dropped the subscript "0" so that (f>o(r,9) — > <j)(r,9)) 

V 2 y = 8ttto 2 (W>* , (1.4) 

--^-V 2 <t> + mV4> = E<f>, (1.5) 
2to 

where V is the Newtonian gravitational potential, and E is an energy eigenvalue. The scalar field 
cf>o(r,9) is then expanded in associated Legendre functions, P[ k (9): 

<j>o{r,e) = 4=Y J Ri{r)P?{9), (1.6) 
^ 4?r i=k 

Similarly, the potential V(r 5 9) is assumed to have no (^-dependence, and thus can be expanded as 

oo 

V(r,9)=Y / V l (r)P l (9). (1.7) 
1=0 

Eqs. H1.6f) and (|1.7fl are then substituted into l|1.4|l and (|1.5fl , and multiplied by P; (9) to obtain a 
system of equations for any particular value of Iq (using the orthogonality of the associated Legendre 
functions). Using this strategy, Silveira and de Sousa were able to obtain solutions for each specific 
combination of Iq and k. A main difference of these solutions relative to the spherically symmetric 

3 Note that 11.31 is clearly not the most general ansatz we could make for stationary solutions of the Einstcin-Klein- 
Gordon system. For example, we could consider p independent bosonic fields <pi(t,x),i = 1,2, • ■ ■ ,p, each satisfying 
ansatz 11.31 for specific values of (k, ui) = [hi, uii), i.e. <f>i (t,x) = <f>i(r, Q)e l ^ u ' it+ki ^ , with a total stress energy tensor 
given by = J^iLi ^iui with each of the given by 12.441 . The stationary solutions x; ki,uji(kj)) would 
then be labelled by p quantum numbers (fej) with p associated eigenvalues (u)j(fcj)), and the spectrum would still be 
discrete. 

4 By which we mean that both the gravitational and matter fields are treated using non-relativistic equations of 
motion. 
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ones, is that the scalar field vanishes at the origin, and hence the rotating star solutions have 
toroidal level surfaces of the matter field, rather than spheroidal level surfaces as in the spherical 
case. 

In 1996 Schunck & Mielke used the ansatz (|1.3[1 to construct rotating boson star solutions. 
Specifically, they chose some particular values of k, specifically k = 0, 1 • • • , 10 and k = 500 and 
showed that solutions to the fully general relativistic equations did exist for those cases. They also 
showed that the angular momentum, J, and the total particle number, N, of the stars are related 
by 

J = kN, (1.8) 

where k is the integer defined in (|1.3[) . An important implication of the above equation is that if we 
consider equilibrium configurations with the same total particle number N, then the total angular 
momentum has to be quantized. This property is in clear contrast with that of a perfect fluid star. 

Later, in the work most relevant to the study of rotating boson stars described in this thesis, 
Yoshida & Eriguchi used a self-consistent-field method j2H] to obtain the whole family of 
solutions for k = 1, as well as part of the family for k = 2. The maximum mass they found for the 
k = 1 case was 1.314Mpj/m. However, their code broke down before they could compute the star 
with maximum mass for the k — 2 case. 

In Table H~l. 31 we summarize some similarities and differences between rotating fermion stars 
and rotating boson stars. For further background information on boson stars we suggest that 
readers consult the review by Jetzer 0], or the more up-to-date survey by Schunck & Mielke [TU|. 



Relativistic Rotating Fermion Stars Relativistic Rotating Boson Stars 

Similarities 

Come in families of solutions parametrized by a single value 
Each family has a maximum possible mass 
Dijferences 

Parametrized by p(0) Parametrized by Ic^'^O)! 

Spheroidal level surfaces of rotating matter field Toroidal level surfaces of matter field for k = 1, 2, • ■ • 
Finite size, abrupt change in d r p, d r p at surface Exponential decay to infinity 

Angular momentum can vary continuously Angular momentum quantized: k = 1, 2, • ■ • 

Table 1.1: Similarities and differences between relativistic rotating fermion stars and relativistic 
rotating boson stars. 



1.2 Critical Phenomena in Gravitational Collapse 

Over the past decade, intricate and unexpected phenomena related to black holes have been dis- 
covered through the detailed numerical study of various models for gravitational collapse, starting 
with Choptuik's investigation of the spherically symmetric collapse of a massless scalar field |7|. 
These studies generally concern the threshold of black hole formation (a concept described below), 
and the phenomena observed near threshold are collectively called (black hole) critical phenomena, 
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since they share many of the features associated with critical phenomena in statistical mechanical 
systems. The study of critical phenomena continues to be an active area of research in numerical 
relativity, and we refer the interested reader to the recent review article by Gundlach for full 
details on the subject. Here we will simply summarize some key points that are most germane to 
the work in this thesis. 

To understand black hole critical phenomena, one must understand the notion of the "threshold 
of black hole formation" . The basic idea is to consider families of solutions of the coupled dynamical 
equations for the gravitational field and the matter field that is undergoing collapse (the complex 
scalar field, </>, in our case). Since we are considering a dynamical problem, and since we assume 
that the overall dynamics is uniquely determined by the initial conditions, we can view the families 
as being parametrized by the initial conditions — variations in one or more of the parameters that 
fix the initial values will then generate various solution families. We also emphasize that we are 
considering collapse problems. This means that we will generically be studying the dynamics of 
systems that have length scales comparable to their Schwarzschild radii, for at least some period 
of time during the dynamical evolution. We also note that we will often take advantage of the 
complete freedom we have as numerical experimentalists to choose initial conditions that lead to 
collapse, but which may be highly unlikely to occur in an astrophysical setting. 

We now focus attention on single parameter families of data, so that the specification of the 
initial data is fixed up to the value of the family parameter, p. We will generally view p as a 
non-linear control parameter that will be used to govern how strong the gravitational field becomes 
in the subsequent evolution of the initial data, and in particular, whether a black hole forms or not. 
Specifically, we will always demand that any one-parameter family of solutions has the following 
properties: 

1. For sufficiently small values of p the dynamics remains regular for all time, and no black hole 
forms. 

2. For sufficiently large values of p, complete gravitational collapse sets in at some point during 
the dynamical development of the initial data, and a black hole forms. 

From the point of view of simulation, it turns out to be a relatively easy task for many models 
of collapse to construct such families, and then to identify 2 specific parameter values, p~ (p + ) 
which do not (do) lead to black hole formation. Once such a "bracket" has been found, 

it is straightforward in principle to use a technique such as binary search to hone in on a critical 
parameter value, p* , such that all solutions with p < p* (p > p*) do not (do) contain black holes. 
A solution corresponding to p = p* thus sits at the threshold of black hole formation, and is known 
as a critical solution. It should be emphasized that underlying the existence of critical solutions 
are the facts that (1) the end states (infinite-time behaviour) corresponding to properties 1. and 
2. above are distinct (a spacetime containing a black hole vs a spacetime not containing a black 
hole) and (2) the process characterizing the black hole threshold (i.e. gravitational collapse) is 
unstable. We also note that we will term evolutions with p < p* subcritical, while those with p > p* 
will be called supercritical. 

Having discussed the basic concepts underlying black hole critical phenomena, we now briefly 
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describe the features of critical collapse that are most relevant to the work in this thesis. 

First, critical solutions do exist for all matter models that have been studied to date, and for 
any given matter model, almost certainly constitute discrete sets. In fact, for some models, there 
may be only one critical solution, and we therefore have a form of universality. 

Second, critical solutions tend to have additional symmetry beyond that which has been adopted 
in the specification of the model (e.g. we will impose spherical and axial symmetry in our calcula- 
tions). 

Third, the critical solutions known thus far, and the black hole thresholds associated with them, 
come in two broad classes. The first, dubbed Type I, is characterized by static or periodic critical 
solutions (i.e. the additional symmetry is a continuous or discrete time-translational symmetry), 
and by the fact that the black hole mass just above threshold is finite (i.e. so that there is a 
minimum black hole mass that can be formed from the collapse). The second class, called Type 
II, is characterized by continuously or discretely self-similar critical solutions (i.e. the additional 
symmetry is a continuous or discrete scaling symmetry), and by the fact that the black hole mass 
just above threshold is infinitesimal (i.e. so that there is no minimum for the black hole mass that 
can be formed). The nomenclature Type I and Type II is by analogy with first and second order 
phase transitions in statistical mechanics, and where the black hole mass is viewed as an order 
parameter. 

Fourth, solutions close to criticality exhibit various scaling laws. For example, in the case of 
Type I collapse, where the critical solution is an unstable, time-independent (or periodic) compact 
object, the amount of time, r, that the dynamically evolved configuration is well approximated by 
the critical solution per se satisfies a scaling law of the form 

t(p) j\n\p~p*\, (1.9) 

where 7 is a universal exponent in the sense of not depending on which particular family of initial 
data is used to generate the critical solution, and ~ indicates that the relation l|1.9|) is expected to 
hold in the limit p — > p* . 

Fifth, and finally, much insight into critical phenomena comes from the observation that al- 
though unstable, critical solutions tend to be minimally unstable, in the sense that they tend to 
have only a few, and perhaps only one, unstable modes in perturbation theory. In fact, if one 
assumes that a Type I solution, for example, has only a single unstable mode, then the growth 
factor (Lyapunov exponent) associated with that mode can be immediately related to the scaling 
exponent 7 defined by (|1.9(l . 

In this thesis we will be exclusively concerned with Type I critical phenomena, where the 
threshold solutions will generally turn out to be unstable boson stars. Previous work relevant to 
ours includes studies by (1) Hawley and Hawley & Choptuik of boson stars in spherically 
symmetry, (2) Noble [3U| of fluid stars in spherical symmetry and (3) Rousseau [31] of axisymmctric 
boson stars within the context of the conformally flat approximation to general relativity. Evidence 
for Type I transitions have been found in all three cases. 
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1.3 Layout 

The remaining chapters of this thesis are organized as follows. In Chap. [21 we summarize the 
mathematical formalism used in the work of this thesis. This includes a brief summary of the 
mathematical model of spacetime, in which the key ingredient to be used is the Einstein field 
equation. We then summarize the ADM (3+1) formalism, which will be used in the study of 
boson stars in spherical symmetry, as well as the (2+l)+l formalism, which is used in the study 
of boson stars in axisymmetry. We also describe the Einstein-Klein-Gordon system, which is the 
fundamental set of PDEs underlying all of our studies. 

In Chap. we summarize the numerical methods used in the thesis. This includes finite dif- 
ferencing techniques which are central to all the calculations shown; the multigrid method, which 
is used in the construction of rotating boson stars, as well as in the solution of the constraint 
equations in the dynamical study of boson stars in axisymmetry; adaptive mesh refinement, which 
is essential in the study of critical phenomena in axisymmetry; excision techniques which are used 
in the study of boson stars in spherical symmetry; and the technique of spatial compactification, 
which is used in the construction of rotating boson stars. 

In Chap. 0] we study Type I critical phenomena of boson stars in spherical symmetry. This 
research can be viewed as an extension of the work reported by Hawley & Choptuik in |5| . Our 
principal new result is compelling numerical evidence for the existence of oscillatory final states of 
subcritical evolutions. We also perform perturbation analyses and show that the simulation results 
agree very well with those obtained from perturbation theory. We then present a rudimentary, but 
stable and convergent implementation of the black hole excision method for the model. Supercritical 
simulations using excision show that the spacetimes approach static black holes at late time, so 
there is no impediment to very long run times (in physical time). 

In Chap. [5] we describe a study of boson stars in axisymmetry. We first present an algorithm 
to construct the equilibrium configurations of rotating boson stars that is based on the multigrid 
technique. We argue that our method is more computationally efficient than methods previously 
used and reported in the literature. More importantly, we obtain numerical solutions in the highly 
relativistic regime for an angular momentum parameter k — 2. We then discuss studies of the 
dynamics of boson stars in axisymmetry. Following Choi's |32| work in the Newtonian limit, we 
show that solitonic behaviour occurs in the head-on collision of boson stars with sufficiently large 
relative initial velocities. We also present a study of Type I critical phenomena in the model. The 
two classes of simulations (collisions of boson stars, and boson stars gravitationally perturbed by a 
massless real scalar field) provide evidence for Type I black hole transitions, as well as scaling laws 
for the lifetime of near-critical configurations of the form H1.9(l . This represents the first time that 
Type I behaviour has been observed in the context of axisymmetric collapse. 

The final chapter summarizes the results, gives overall conclusions and points to some directions 
of future work. Several appendices providing various technical details are also included. 
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1.4 Conventions, Notation and Units 

Throughout this thesis the signature of metric is taken to be ( h + +). Spacctimc indices of four 

dimensional (1 temporal + 3 spatial) tensors are labeled by lower case Greek letters (a, (3, 7, • • • ). 
Spatial indices of three dimensional (3 spatial) tensors are labeled by lower case Latin letters starting 
from i (i,j,k,---). Spacetime indices of three dimensional (1 temporal + 2 spatial) tensors are 
labeled by the first few Latin indices (a, b, c, • • • , h), and spatial indices of two dimensional (2 
spatial) tensors are labelled by upper case Latin letter (A,B,C,---). The Einstein summation 
convention is implied for all types of indices. That is, repeated (one upper and one lower) indices 
are automatically summed over the appropriate range. Covariant derivatives are denoted by V or 
by a semi-colon ";" • Ordinary partial derivatives are denoted by d or by a comma "," • Conventions 
for the Riemann and Ricci tensors are V^V^jV^ = \Rafj-ysV 5 , R^ v — R a ^ a u- 

We adopt a system of "natural" units in which G = c = K = 1. In this system the unit time, unit 
length and unit mass are known as the Planck time, Planck length and Planck mass, respectively. 
Specifically, we have 



T p i = W-g- = 5.39 x l(T 44 s, (1.10) 



L p i = = 1.62 x ltr 35 m, (1.11) 

l~hc 

M p i = J— = 2.18 x l(T 8 kg. (1.12) 
V (j 

The corresponding energy scale is M v \ c 2 w 10 19 GeV. For the scalar field, the action has dimension 
[h]. Hence [d 4 x(d(/)) 2 ] = [h]. Therefore, [<j>\ = [Vh]/L = y/M/T. Thus, <p is measured in units of 
VM P i/T pl = 6.36 x 10 17 v/kg7I 

Finally, we adopt a nomenclature commonly used in numerical work, whereby we refer to 
calculations requiring the solution of partial differential equations in n spatial dimensions as "nD" . 
In particular, we will often refer to spherically symmetric computations (whether time-dependent 
or not) as "ID" , and axisymmetric calculations (again, irrespective of any time-dependence) as 
"2D". 
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CHAPTER 2 

MATHEMATICAL FORMALISM 

In this chapter we briefly summarize the mathematical description of spacetime, as well as the 
ADM (or 3+1) and (2+l)+l decompositions of spacetime. These decompositions — which form the 
bases for the numerical calculations described in the remainder of this thesis — are well described in 
the literature, and we therefore restrict our discussion here to a summary of the physical picture, 
and statements of the key equations that will be used in subsequent chapters. 

2.1 Mathematical Model of Spacetime 

In general relativity the elementary entities are events, each of which is characterized by a time, 
i, and location, x l , at which the particular event occurs. According to Einstein's geometric view 
of the gravitational interaction, gravitational physics is concerned with the relationship between 
these events: in other words, it describes the "structure" of a collection of events M. = {(t,x 1 )}. 
Specifically, we are to view the collection of events as a surface with a certain structure imposed on 
it. In mathematical language, the collection of events forms a manifold, where manifold is simply a 
generalization of the concept of surface, and the structure that describes the gravitational physics 
turns out to be a metric on the manifold. The pair of objects, the manifold, A4, and the metric, 
g llVl comprise spacetime, which is of primary interest in general relativity. The study of spacetime 
thus becomes the study of the explicit form of the metric, or, equivalently, of the geometry of 
spacetime. 

The physics that we are looking for should also tell us how the distribution of matter affects the 
structure of spacetime. In other words, the theory should have a way of prescribing matter fields, 
and the relationship between these matter fields and the geometry of spacetime. Mathematically, 
we describe the matter as some tensor fields on the manifold, and the relationship between the 
matter and the geometry is written in the form of field equations. We will make some further 
assumptions about the spacetime so that some general physical principles — of whose validity we 
are confident — will be satisfied. Following Hawking and Ellis [HSj; these assumptions are: 

1. Local causality: special relativity posits that no signal can travel faster than the speed of light. 
Therefore, equations for any matter fields should respect this property. Causality effectively 
partitions the neighborhood of any point p into two sets: those points which are causally 
connected to p and those which are not. This assumption allow us to determine the metric 
up to a conformal factor (i.e. up to an overall location-dependent scale). 

2. Local conservation of energy and momentum: conservation laws generally reflect symmetries 
in physical laws. Energy and momentum conservation — related to invariance of physical laws 
under temporal and spatial translations — is a cornerstone of physics, with an extremely solid 
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empirical basis. Conservation of energy and momentum can be expressed as a condition 
on the energy- momentum tensor, T),„, which describes the matter content of the spacetime. 
Knowledge of this tensor determines the conformal factor. In general the conservation equa- 
tions would not be satisfied for a connection derived from a metric g = fl 2 g. Observing two 
paths j(t) and j'(t) of small test particles would allow us to determine O up to a constant 
factor, which simply reflects one's units of measurement. 

3. Field equations: this is the most stringent assumption that we make, but it also yields the 
key equation that we wish to solve, namely the Einstein field equation: 



Here G^ Vl and R are, respectively, the Einstein tensor, the Ricci tensor and the Ricci 
scalar. (We note that, in contrast to [33| . we further assume that the cosmological constant 
vanishes.) We observe that predictions made on the basis of this equation agree, within 
experimental/observational uncertainties, with all experiments and observations that have 
been performed to date. (A survey of experimental tests of general relativity can be found 
in |34| . and an up-to-date review of the state of the art in experimental verification of the 
theory is given in |35l K-t6|). The Einstein equation is the mathematical expression of the 
statement that "matter tells spacetime how to curve" . In turn, the curvature of spacetime 
"tells matter how to move" through the appearance of the metric, and gradients of the metric 
in the equations of motion for the matter fields. 1 



One of the great conceptual revolutions of relativity is the unification of the concepts of space 
and time into a single entity. This unification is partly motivated by the fact that, under general 
coordinate transformations, time and space "mix" together, and it is thus no longer meaningful to 
talk about absolute space or absolute time. In particular, the field equation (|2.1J) . which is invariant 
under general coordinate transformations, relates the entire spacetime geometry to the distribution 
of matter-energy in space and time. The value of physics, on the other hand, is often rooted in its 
ability to predict what happens in the future, given certain information that is known at some initial 
time (where time must generally be defined with respect to a given family of observers). In other 
words, one traditional way of solving a physics problem is through the study of the dynamics of the 
system under consideration; i.e. given the state of a physical system at a specific instant of time, 
we use dynamical equations of motion to evolve the state to yield information about the system 
at future (or past) times. This approach of viewing the solution of (|2.1|l as the time development 
of some initial data is known as the Cauchy problem for relativity. Note that realization of this 

1 We note that Einstein gravity can be formulated in a coordinate-free (i.e. geometrical fashion); this in fact is 
a necessary consequence of the general covariance of the theory. Nonetheless to construct spacetimes (i.e. to solve 
Einstein's equations, particularly via numerical means, we must adopt specific coordinate systems (also known as 
charts or an atlas). Although this can be a subtle issue, we will assume in the following that the coordinate systems 
we choose cover the region of spacetime of interest (i.e. we adopt global coordinate systems). 




(2.1) 



2.2 The ADM (3+1) Formalism 
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strategy will in some sense undo Einstein's feat of unification, i.e. it will split up spacetime into 3 
"space" dimensions plus 1 "time" dimension — hence the terminology "3+1" split. The split will 
involve the introduction of some specific coordinate system, including a definition of time, and a 
prescription of initial data at the specified initial time. The initial data are to be chosen so that, 
as much as possible, the resulting evolution will accurately model the physical system of interest. 
We will then march the data forward in time using the dynamical form of the field equation. Our 
task then, is to rewrite equation (|2.1|) as a set of partial differential equations (PDEs) suitable for 
solution as a dynamical problem. As it turns out, the PDEs themselves naturally decompose into 
two classes: one which involves only quantities defined at a given instant of time (the constraints), 
and another that describe how the dynamical variables evolve in time (evolution equations). 

The split of spacetime into space and time as sketched above, was first introduced by Arnowitt, 
Deser and Misner |371 13*5| (hence the name ADM formalism) . The original motivation for the 
development of the formalism was, in fact, the desire to cast the field equations in a form suitable 
for quantization. Although this so-called canonical quantization of the gravitational field remains 
a largely unsolved problem, the ADM approach has proven to be a very useful basis for classical 
computations in numerical relativity. 

The explicit 3+1 decomposition proceeds as follows. We suppose that spacetime is orientable, 
which, roughly speaking means that we can time order the events in spacetime. We then slice 
up the spacetime into different slices of constant time, and label each slice with a value of the 
time parameter, t. We demand that the slicing is such that any two events on any given slice are 
spacelike-separated; i.e. that the slices are spacelike hyper surfaces 2 . The collection of these slices 
{£ t } becomes a particular "foliation" of spacetime. Mathematically this foliation is described, at 
least locally, by a closed one form, fl — dt (the closed nature of the form enables us to use f as a 
unique label for the slice). The norm (length) of this form quantifies how far (proper distance) it 
is from one slice, E 4 , to a nearby slice, Y, t +dt- 

||0|| 2 = g^Vl^ u = -oT 2 . (2.2) 

Here, the positive function, a, is called the lapse function and encapsulates one of the four degrees 
of coordinate freedom in the theory. The unit-normalized one form, 

n M = -af2 M , (2.3) 

and the associated unit-norm vector, — g^ v n v , then allow us to define a projection operator, 

J." v = 5" v + n»n v . (2.4) 

Note that since n M is timelike, we have g fl ^n ti n 1 ' — g^n^n^ = n^n^ = — 1. Using the projection 
operator, we can construct tensors that "live in" the hypersurfaces by projecting all of the indices of 
generic 4-dimensional (spacetime) tensors. For example, for a rank-3 tensor with one contravariant 

2 Roughly speaking, a hypersurface is a d — 1 surface in a (i-dimensional manifold. More precisely, it is the image 
of an imbedding of a d — 1 dimensional manifold. 
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and two covariant indices, we can define 

j_ T a M =± a K ±% \y 7 rv . (2.5) 

By construction, any tensor denned in this fashion is orthogonal to the unit normal, so that, for 
example, we have n a _L T Q ^ 7 = -nP _L T a ^ 7 = n 1 _L T a j a 7 = 0. Hence such a tensor naturally lives 
in the 3-dimensional hypersurfaces, and is known as a spatial tensor. In particular, we construct 
the induced (or spatial) metric, 7^, via 

7 M „ =_L g^u = g^u + n^n u , (2.6) 

and note that 7^ describes the intrinsic geometry of the hypersurfaces. Similarly, we can use 
projection to define a 3-covariant derivative, D ai that is compatible with the 3- metric, 7^„ (i.e. so 
that -D Q 7 M „ = 0). For example, for a spatial vector, (V^n^ — 0), we have 

D a V» =± V a V» =±P a ±" vVpV" , (2.7) 

where V Q is the spacetime covariant derivative. Using D a and standard formulae of tensor calculus, 
we can then compute the various 3-dimensional curvature tensors (Riemann tensor, Ricci tensor, 
Ricci scalar etc.) that characterize the geometry of the hypersurfaces. 

In the 3+1 approach, the components of the 3- metric, 7^, are to be viewed as dynamical 
variables of the theory; the geometry of spacetime then becomes the time history of the geometry 
of an initial hypersurface. In order to characterize how the geometries of two nearby slices differ, 
we must consider the manner in which the slices are embedded in the enveloping spacetime. This 
information is encoded in the (symmetric) extrinsic curvature tensor, which, up to a sign, is simply 
the projected gradient of the unit normal to the hypersurfaces: 

= K vtl = - ± V^riv . (2.8) 

From the above definition, it is clear that K^ v is a spatial tensor, and thus can also be viewed as 
living in the hypersurfaces. Roughly speaking, the components of the extrinsic curvature can be 
thought of as the velocities (time derivatives) of the 3- metric components (see equation (|2.17|0 . 
so that 7^1/ and K^ v are, loosely, dynamical conjugates. With the above definitions of the spa- 
tial metric and the extrinsic curvature (also known as the first and second fundamental forms, 
respectively), one can derive the Gauss-Codazzi equations 3 

-L RaPfS = R a 0<y5 + KspK-ya — KffiKsa , (2-9) 

-L RuPin = DpK a7 — D a K/3 7 , (2.10) 

where ^Rafi^s denotes the Riemann tensor constructed from the spatial metric, 7^, and where, 
following York |38j . a n in a covariant index position denotes that the index has been contracted 
with n M — for instance, Wn = W^n* 1 . (However, and again following York, a n in a contravariant 

3 See, for instance, 1331 . 
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index position denotes contraction with — n^.) 

The Einstein equation, together with the following definitions: 

p = T flh = T^n\ (2.11) 
f = ir" = -l {T^nv) , (2.12) 

where T^ v is again the stress energy tensor, allow us to rewrite the Gauss-Codazzi equations as 

(3) i? + K 2 - K i:j K ij = 167rp , (2.13) 



DjK ij - D l K = 8nf . (2.14) 

Here, is the 3-dimensional Ricci scalar, K = 7 y K+j is the trace of and we remind the 
reader that Latin indices range over the spatial values 1, 2, 3. We also note that spatial indices on 
spatial tensors are raised and lowered using ^ and 7y, respectively, and that 7y7 J = 5 l k- 

(|2.13l) and (|2.14l) involve only spatial tensors and spatial derivatives of such tensors, and hence, 
from the point of view of dynamics, represent constraints that must be satisfied on each slice, 
including the initial slice. They are commonly called the Hamiltonian constraint and momentum 
constraint respectively, and are a direct consequence of the four-fold coordinate freedom of the 
theory. 

So far we have used only one coordinate degree of freedom, which was in our specification 
of the manner in which we foliated the spacetime into a family of spacelike hypersurfaces — or, 
equivalently, in the specification of the lapse function, a, at all events of the spacetime. For any 
fixed time slicing, however, we are allowed to relabel the spatial coordinates of events intrinsic to 
any of the slices as illustrated in Fig. 12.11 As shown in that figure, it is convenient to describe 
the general labeling of events with spatial coordinates relative to the specific labelling one would 
get by propagating the spatial coordinates along the unit normal field, n". In general, the spatial 
coordinates, will be "shifted" relative to normal propagation, and this shift will be quantified by a 
spatial vector, known, naturally enough, as the shift vector. An observer who is at rest with 
respect to the coordinate system illustrated in Fig. 12.11 will thus move along a worldline with a 
tangent vector given by 

f = an" + 0" , (2.15) 

where we again emphasize that P^n^ = 0. 

The specification of the three independent components of the shift vector, along with the lapse 
function exhausts the coordinate freedom of general relativity, and defines a coordinate system 
(t, a; 4 ) (sometimes called a 3+1 coordinate system) that we will assume covers the region of space- 
time in which we are interested. In such a coordinate system, and using a four-dimensional gener- 
alization of the Pythagorean theorem, we can then write the spacetime displacement in 3+1 form 



ds 2 = {-a 2 + f3 t p l ) dt 2 + 2(3, dt dx l + 7li dx l dx j 



(2.16) 
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Figure 2.1: This figure illustrates a schematic representation of the ADM, or 3+1, decomposition. 

a is the lapse function and encodes the proper distance, adt, between two infinitesi- 
mally separated hypersurfaces, E(£) and £(£ + d£). j3 % is the shift vector and describes 
the shifting of spatial coordinates, x l , relative to normal propagation. PP" is the 
vector £ M = an^ + (3^, where n M is the unit timelike normal to the hypersurfaces, and 
is proportional to the 4-velocity of an observer with fixed spatial coordinates. 

From the definitions of the vector field, £ p Q2.15[l . and the extrinsic curvature (|2.8(l . one can show 
that 

£ tll3 = -2aKij + Dify + Djfa , (2.17) 

where £ t denotes Lie differentiation along £ M (in a 3+1 coordinate system £ t reduces to partial 
differentiation <9 t ). This last expression can be viewed as a set of evolution equations for the 3-metric 
components, and makes the interpretation of the extrinsic curvature components as "velocities" of 
the 7ij manifest. 

Making the further definitions 

V = -LI^, (2.18) 
S = fJS l3 , (2.19) 

the remaining components of the Einstein equation can be written as 

£ t K i:j = - DiDj a + a | (3) Rij - 2K lk K k 2 + K^K - 8tt / 
+p k D k K, tJ + K lk D^ k + K kj DiP k , 

where ^ Rij is the 3-dimensional Ricci tensor. These last equations are to be viewed as evolution 
equations for the extrinsic curvature components. 

In summary, the 3+1 decomposition applied to the Einstein field equation l|2.1(l yields a set of 



(2.20) 
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4 constraint equations, (|2.13() . I|2.14[l . and 12 (first order in time) evolution equations, (|2.17|l and 
(|2.2U|) . We note that, as a consequence of the contracted Bianchi identity (which itself follows from 
the coordinate invariance of the theory), the evolution equations guarantee that data that satisfies 
the constraints at some instant in time (at the initial time, for example), is evolved to data that 
satisfies the constraints at any future or past times. 

We conclude this section by observing that a key idea in the 3+1 approach is to use a timclikc 
vector field, (whose existence is guaranteed by the Lorentzian signature of the spacetime metric), 
to construct a projection operator that is then used to define quantities on surfaces orthogonal to 
n M (surfaces of simultaneity), and to derive equations governing those quantities from the covariant 
field equations. This same idea is used in the (2+l)+l formalism described in the next section, 
the only difference being that we will now decompose with respect to a spacelike vector field, £ M , 
which is a Killing vector field. 

2.3 The (2+l)+l Formalism 

Suppose we stack a collection of identical planar figures — such as the text string "(2+l)+l" — to 
create a 3-dimensional structure as shown in Fig. 12.21 For the purposes of discussion, we will think 
of the 3-dimensional object as being of infinite extent in the z-direction. To describe this structure 
it suffices to state that its z — const, cross section is the planar text "(2+l)+l", and that it has a 
symmetry in the z direction. Similarly if there is a symmetry in spacetime, we can state what the 
symmetry is, then focus on the detailed structure of a "cross section" orthogonal to the symmetry 
direction. This technique of "dividing out" the symmetry and then restricting attention to the 
cross section, or projected space, (or, in mathematical parlance, the quotient space 4 ), underlies the 
(2+l)+l formalism. In particular, if the spacetime in which we are interested is axisymmetric, we 
can first project quantities and equations onto a "plane" perpendicular to this symmetry, and then 
perform the equivalent of an ADM space-plus-time decomposition on the dimensionally reduced 
spacetime. 

The (2+l) + l formalism described here is based on the work of Maeda and his collaborators |39j . 
and, as just stated, combines (1) the method of dimensional reduction (dividing out the symmetry, 
or projection along the symmetry) of a spacetime possessing a Killing vector field, as originally 
developed by Geroch |40|. and (2) a 3-dimensional analogue of the ADM approach introduced in 
the previous section (Sec. 12. 2|) . Other references for the application of this formalism to problems 
in numerical relativity include |41l |H] . 

As previously mentioned, we are interested in the case of axisymmctry, and thus assume that 
the spacetime under consideration contains a spacelike Killing vector field, with closed orbits. 
Choosing an azimuthal coordinate, ip, adapted to the symmetry, the Killing vector can be written 

4 The equivalent relation can be denned as p ~ q iff p and q lie on the same integral curve of the Killing vector 
field £. 
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Figure 2.2: A 3-dimensional object with continuous symmetry, which can be described as having 
a cross section that is the planar text string "(2+l)+l", and which has a symmetry 
along the z direction. 



as: 



(2.21) 



To perform the projection, we first construct a projection operator, H a p, similar to that previ- 
ously used in the 3+1 decomposition 5 : 



H a p = 5° 



(2.22) 



Using this operator, we can construct tensors intrinsic to the cross sections via projections of 
spacetime tensors. In particular, we can project the spacetime metric itself (or equivalently, lower 
the contravariant index of H a p) to yield the metric, H a p, on the dimensionally reduced space: 



~H a p = g a p 



(2.23) 



Again paralleling the 3+1 development, we can also use projection of the spacetime covariant 
derivative, V a , to define an induced covariant derivative operator, D a , in the cross sections that 
satisfies D a 1-t^ v = 0. 

Since it is a (Lorcntzian) metric on a 3-dimensional manifold, TL a p clearly has fewer degrees of 
freedom than the spacetime metric, g a p. The "missing" degrees of freedom can be conveniently 
encoded in the norm, s, and the twist, u> a , of the Killing vector, £ M : 

5 This definition of projection operator, which is constructed from a vector (£) instead of a 1-form (Q) makes a 
significant difference for the relation between 4-dimensional objects and 3-dimensional ones. It is straightforward to 
get the contravariant objects in the former situation, while it is straightforward to get the covariant objects in the 
latter case. The main reason is that £ has only 1 non-zero component in the contravariant form, but f! has only 1 
non-zero component in the covariant form. For instance, ( 3 )y a ~ -LjJ ( 4 )yM = (5« — £ M £ a /s 2 )( 4 - > V' M = ( 4 )y a , since 
C a =0. 
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,s 2 = (2.24) 
io a = e a PlS e^- (2-25) 

where e a p~ f $ is the usual Levi-Civita symbol. Note that since w a £; Q = L^oj a = C^s 2 = 0, s 2 and u> a 
are tensors intrinsic to the cross sections. For notational convenience, we further define 



,2 



w a u) a , (2.26) 



A = s 2 . (2.27) 

The basic equations for a spacetime with a Killing vector field can now be written as 6 : 

D [a tu fJ] = -e a ^ s CR S e^ e , (2.28) 
D a u a = ^ D^X, (2.29) 

D 2 \ = ^D a XD a X-^-2R al3 Cf , (2-30) 

n a() = (uj a ui - n a0 uj 2 ) + - s D a D ()S + n» a n v (3 R^ , (2.31) 

where H a p denotes the Ricci tensor associated with the induced metric Tt a f3- Note that H2.28(l - 
(|2.31|1 completely describe the dimensional reduction piece of the (2+l)+l formalism, and are 
written in purely geometrical terms. To incorporate matter fields, we can manipulate (|2.28(l - (|2.31() 
and use the Einstein equation to get: 



D[ a uj b ] = 8irse abc T c , (2.32) 



UJa 
*3 



= 0, (2.33) 



D 2 s = [tg-T , (2.34) 



-D 2 s = -4-4.^ 
s 2s 4 V s 2 

Gab = Ti-ab — ^Hab^ = ^T ab , (2.35) 

where 1Z = H ab lZ a b- Here we note that indices a, 6, • • • range over the temporal and (two) spatial 

6 See appendix of 1401 for a detailed derivation. 
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dimensions of the reduced space, and that we have defined the following quantities: 



Tab = rV a ri 5 bTj5 , 



ni C rpQL 

n a 1 q 



T — 1 U ( T <f"t 
■lab — Tab — ^rtab \T -J" 



J_ j U a U) b 1 

8tt I 2s 4 s 



D a D b s - -H ab D's 



(2.36) 
(2.37) 
(2.38) 

(2.39) 



Note that the 4-dimensional Einstein equations have been effectively replaced by l|2.32|l - H2.35|l . 
where the first two equations determine the twist ui a 7 , the third equation determines the norm, s, 
of the Killing vector field, and the fourth governs the reduced metric, TL a p. 

The last stage in the derivation of the (2+l) + l equations involves a space-plus-time split 
(completely analogous to the 3+1 split) of the reduced field equations. Again, the interested 
reader is referred to [HHI f° r full details. 



2.4 Boson Stars — the Einstein-Klein-Gordon System 

In this section we briefly describe the theoretical framework for boson stars, which are equilibrium 
configurations of a massive, complex Klein-Gordon field coupled to gravity — we will thus refer to 
the model as the Einstein-Klein-Gordon system. 

We derive the equations for the Einstein-Klein-Gordon system using a variational approach. 
The model is described by an action 

S = J d 4 x^(C g + C^) , (2.40) 



and 



where 

C = ±R, (2-41) 

= ~ (V^V M 0* + m 2 |0| 2 ) . (2.42) 

Here R is the spacetime Ricci scalar and <p is a massive complex scalar field with bare (particle) 
mass to. Variation of the action with respect to the metric g^v yields the Einstein equation 

~ \g^R = 87rT M , , (2.43) 

where 

Tpv = T*, = \ [(V^V„r + V^V^*) - 9liv ( V Q 0V Q r + m 2 |0| 2 )] . (2.44) 

7 Helmholtz's decomposition theorem tells us that a vector field can be determined by its curl and divergence. 
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Similarly, variation with respect to the field <fi gives the Klein- Gordon equation 

V^V M - m 2 (j) = . (2.45) 

The coupled equations (|2.43|) and l|2.45[l constitute the equations of motion for the Einstein-Klein- 
Gordon system, and, in particular, are the fundamental equations used in our subsequent study 
of boson stars. When we (gravitationally) perturb boson stars using an additional real, massless 
scalar field, 03, we simply add an extra term C$ 3 to the total Lagrangian density 

C = C g +£^+C 4>3 , (2.46) 

with 

£03 = -^V M 03. (2.47) 
This yields an additional contribution, Tj%j>, to the total stress energy tensor 

= T* v + T# , (2.48) 

where 

?il = V^V^a ~ V Q 3 V Q 3 , (2.49) 
as well as a field equation (the massless Klein- Gordon equation), for 03 

V^V M 3 = O. (2.50) 



In the following chapters we will consider more explicit forms of the above equations of motion, 
once we have imposed symmetries and fixed specific coordinate systems. 
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CHAPTER 3 



NUMERICAL METHODS 



In this chapter we briefly summarize the numerical methods used in this thesis, as well as some 
continuum strategies central to our numerical calculations. As with the previous chapter, since 
most of the topics discussed here are well documented in the literature, we refer the interested 
reader to appropriate references for additional details. We also note that we use the common 
subscript notation for partial differentiation, e.g. u x = du/dx, in this chapter. 

3.1 Finite Difference Techniques 

The basic technique that we use in our computational solution of model physical problems is 
finite differencing. Finite differencing is a very commonly used method in the numerical solution 
of partial differential equations (PDEs) , and continues to dominate work in numerical relativity. 
The basic idea of the approach is to replace partial derivatives by suitable algebraic difference 
quotients, and the fundamental assumption underlying the efficacy of the technique is that the 
functions to be approximated are smooth so that they can be Taylor-expanded about any point in 
the computational domain. The first step in the development of a finite difference approximation 
(FDA) to a PDE (or more generally, to a system of PDEs) is the introduction of a discrete grid (or 
lattice, or mesh) of points, which replaces the continuum physical domain. We then approximate 
continuum field quantities u = (u,v,w, ■ ■ ■) by a set of grid functions u h = (u h , v , w h , • • • ) which 
will constitute a solution of our FDA. For simplicity of presentation, we will now assume that our 
vector of fields, u, has a single component, u(t, x), which is a function of time and one spatial 
dimension, but the approach and techniques discussed in the following are easily generalized to the 
cases of multiple fields and/or multiple spatial dimensions. 

We will generally restrict attention to so-called uniform grids wherein the spacing between 
adjacent grid points is constant in each of the coordinate directions. Denoting these regular mesh 
spacings in space and time as h and k, respectively, and assuming that we are solving our PDE 
on the continuum domain, x m ; n < x < x max , t > 0, our finite difference grid is given by (xj,t n ), 
where 



xj = x min + (j -l)h, j = !,-■■ , 



N 3 



X ) 



(3.1) 



with 



h = 



x 



max 



- X 



mm 



(3.2) 



N x -1 



and 



t n =nk, n = 0, . 



(3.3) 



For any grid function u h 



€ u 



h , the value at (xj,t n ) is denoted by uj and is an approximation of the 
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continuum value u(xj, t n ). Convergence of the finite difference approximation is then the statement 
that lim^^o.fe^o u j = u{xj,t n ). We remark at this point that, operationally, investigation of the 
behaviour of a finite difference solution, ii 11 , as a function of the mesh spacings h, k (holding other 
problem parameters fixed) provides a very powerful and general methodology for assessing the 
accuracy of the solution. We have used this strategy of convergence testing extensively in assessing 
the correctness and accuracy of the results described in later chapters. 1 In addition, for all of our 
time-dependent calculations, when we vary the spatial mesh spacing(s) during a convergence test, 
we also vary the time step, k, so that (for the one dimensional case), the ratio A = k/h is held 
constant (A is often called the Courant factor, or Courant number). Thus, our FDAs tend to be 
characterized by a single overall discretization scale. 

We now sketch one general technique that can be used to derive finite difference formulae. As 
mentioned above, we want to use difference quotients (algebraic combinations of grid-function val- 
ues) to approximate derivatives. For instance, suppose we wish to approximate the first derivative 
u x (x) of our function u(x). Suppressing the time index, n, we Taylor series expand about x = Xj 

u j+ih = [u]j + ih[u x ]j + -(ih) 2 [u xx ]j + 0(h 3 ) , (3.4) 

for i = —I, • • • , —2, —I, 0, 1, 2, • • • , /. Now suppose we want to use a three-point formula (or a 
three-point stencil) to calculate our approximation of the first derivative (see Fig. I3.ljl . 



h h 
• • • 

J- 1 J J+ J 

Figure 3.1: A 3-point finite difference stencil suitable for computation of a second-order approx- 
imation to a first derivative 



We consider a linear combination, 5Zi=— i CiUj+ih, of the grid function values at the grid points 
Xi-i, Xi, Xi+i, and equate it to the first derivative, u x , evaluated at x = xf 



l 



Solving for the coefficients Cj, we find 



^ CiUj +ih ~ [u x ]j . (3.5) 



c = 0, (3.6) 
Cl = + 2^ 

1 Our determination of rotating boson star data in Chap. |S] is an exception. We typically compute there on a 
N r X Ng = 100 X 15 grid. Calculations with, e.g. N r X Ng = 200 X 30 do not converge. We feel this non-convergence 
is due to the regularity problems that we report in that chapter. 
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and hence our FDA is 

u x (x j ) = Uj+1 - h Uj - 1 +0(h 2 ). (3.7) 

Since the truncation error of this formula is 0(h 2 ), we call the approximation second order. Sim- 
ilarly, we can derive expressions for other derivatives (with varying degrees of accuracy) using 
appropriate stencils. For example, we have |42| 

Ux{Xj) = ^±l_ll +0 (h), (3.8) 

= +0(h), (3.9) 

U xx ( Xj ) = "j+l- 2 ^+"j-l +0(fe 2 ); (31Q) 

= ""^ +2 + + U >- 2 + Q(fr 4 ) . (3.11) 

In discretizing time-dependent PDEs we make exclusive use of Crank-Nicholson schemes, with 
second order spatial differences. For the case of one space dimension (1+1 problems), such a 
scheme uses the stencil shown in Fig. 13.21 



Q q n + 1 




-Q O n 

j j+1 



center-point or scheme : x = x ., t = t 



Figure 3.2: Stencil for an 0(h 2 ) Crank-Nicholson scheme for a PDE in one space dimension and 
time. 



The key idea of a Crank-Nicholson method is to keep the differencing centered in time as well as 
in space. As an example, we consider the simple advection equation 



u t (t,x) = au x . 



(3.12) 



where a is a constant. The Crank-Nicholson scheme for the equation is then 

n+l n 



At 



a 
2 



u — u 
2Ax 



n+l n+l 

2Ax 



(3.13) 



and is 0(Ai 2 , Ax 2 ) accurate as can be verified via Taylor series expansion about the center-point 
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of the stencil, (xj, t n+1 / 2 ) (see Fig. 13.211 . Note that the scheme involves the FDA of the spatial 
derivative applied at both the retarded (n) and advanced (n+ 1) time levels. This leads to coupling 
of the advanced-time unknowns, (i.e. the scheme is implicit), and is crucial for both the 

accuracy and stability of the scheme. 

The Crank-Nicholson approach can be applied in a straightforward manner to any set of time- 
dependent PDEs that have been cast in first-order-in-time form, and for the case of the evolution 
equations treated in this thesis, the discrete equations that result can be solved iteratively in an 
efficient manner. 

For time-dependent FDAs, the notions of convergence and stability are intimately related. We have 
generally not encountered serious stability problems with the difference schemes used below, but 
nonetheless refer the interested reader to 03 for discussion of this important subject. 

3.2 Multigrid Method 

In this section we introduce the basic concepts and techniques of the multigrid method for the 
solution of FDAs arising from the discretization of elliptic boundary value problems. This method 
(including a modification for eigenvalue problems) is used extensively in our construction of station- 
ary solutions describing rotating boson stars, and in the solution of the elliptic constraint equations 
and coordinate conditions that arise in our study of the dynamical evolution of axisymmetric boson 
stars. A more detailed introductory discussion of the method can be found in |441 145| . while |4"o] 
provides an extensive reference source. 

Traditional methods for solving elliptic PDEs, such as successive over-relaxation (SOR), are easy 
to implement, but suffer from slow convergence rates (particularly in the limit of small mesh sizes), 
and have contributed to the general impression that, relative to equations of evolutionary type, the 
finite difference solution of elliptic equations is computationally expensive. However, the multigrid 
method, first introduced by Brandt in the 1970s 07] i is capable of solving discrete versions of 
elliptic PDEs (even systems of nonlinear PDEs) in a time that is proportional to the total number, 
A, of grid points in the discretization — i.e. in O(N) time. The basic idea underlying multigrid is 
to use a hierarchy of grids with different resolutions (mesh spacings) to solve a particular problem 
instead of using a single grid. By an intelligent transfer of the problem back and forth between 
the various grids, one can speed up the convergence rate of the solution process on the finest grid. 
Central to the efficient operation of most multigrid solvers is the observation that straightforward 
relaxation (in particular, Gauss-Seidel relaxation, 45, 48 ) is very good at removing high frequency 
error components, but very bad at removing low frequency ones. At the same time, once the 
solution error has been smoothed by relaxation, we can sensibly pose a coarse grid version of the 
problem at hand, the solution of which can be computed in a fraction of the time required to solve 
the fine grid problem. 

Perhaps the best way of illustrating the multigrid method is by consideration of a simple 
example. Suppose we want to solve the following linear elliptic equation 



Lu = f, 



(3.14) 
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where L is some linear elliptic operator, u is the continuum solution, and / is a source function 
(we assume L is linear strictly to keep the presentation simple; as already mentioned, the multigrid 
technique can also be applied to nonlinear elliptic equations). As in the previous section we replace 
the PDE by a finite difference approximation, defined on a uniform grid characterized by a mesh 
spacing, h 

L h u h = f h . (3.15) 

Here L h is our FDA of L and u h and f h are the grid versions of the solution and source func- 
tion, respectively. Suppose we have an approximate solution (or guess) u h to the above equation. 
Then we consider the difference, v h , between the exact discrete solution, u h , and the (current) 
approximation, u h , to u : 



v h = u h -u h . (3.16) 



We call v h the error or correction, since it is the quantity by which we must correct u , so that 
we get the solution of the finite difference approximation. Solving this last expression for u h , and 
substituting the result in l|3.15|l we have 

L h (u h + v h ) = f h , (3.17) 



L h v h = - (L h u h - f h ) = -r h , (3.18) 

where we have defined the residual or defect, r h . Focus is now shifted to the solution of (|3.18|) for 
the correction, v . Once v h — or more generally, some approximate solution, v , of l|3.18(l — is in 
hand, we can update u h via 

u h :=u h + v h . (3.19) 

Of course, H3.18fl is fundamentally no easier to solve than l|3.15|) . However, if we apply a few 
relaxation sweeps to <|3.18[1 . thus smoothing the residual, r , as well as the (unknown) correction, 
v , we can transfer the problem to a coarser grid with mesh size H, where H is a multiple of h 
(typically H = 2h). That is, we consider the coarse grid problem 

L H V H = _ r H = _ I H r H _ (3 20) 

where the coarse grid source function, r H , is obtained by the application of a suitable restriction 
operator Iff to the fine grid residual, r h (i.e. iff transfers a fine grid function to a coarse grid). 

Particularly for 2D and 3D cases, the coarse grid problem will be significantly less expensive to 
solve than the fine grid problem. Once the solution v H is obtained we can approximate the desired 
correction v h by using an prolongation operator, 1^, which transfers a coarse grid function to a 
fine grid: 

v h = I h H v H , (3.21) 



Once the fine grid function has been updated using the prolonged correction, a few more relaxation 
sweeps are applied to smooth out any high frequency error components that are introduced in the 
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coarse-to-fine transfer of the correction. 

We note that the choice of suitable restriction and prolongation operators will generally depend 
on the discretization that is used, as well as on the specific type of relaxation strategy that is 
employed to smooth the residuals/corrections. More importantly, the above sketched process of a 
two level correction can be easily generalized to a multi-level correction, wherein one uses an entire 
hierarchy of mesh scales (for example h, 2h, Ah, • • • ) and where the cost of solving the discrete 
equations on the coarsest grid can be made computationally negligible. The algorithmic process by 
which the solution is manipulated first on the finest grid, then on the auxiliary coarse grids, and 
finally on the finest grid again is called a multigrid cycle. Although it is sometimes advantageous 
to use other types of cycles, we use V-cycles exclusively in our work, wherein the problem solution 
sequence is If, If — 1, • • • , l c + 1, lc, lc + 1> • • ■ If — 1, If, and where the integers l c and If (If > l c ) 
label the coarsest and finest levels, respectively, of refinement used. 

3.3 Adaptive Mesh Refinement 

In this section we introduce the method of adaptive mesh refinement (AMR). This method is 
important in the finite difference solution of problems that exhibit significant variation in the length 
and time scales that must be resolved (dynamic range), particularly when the specific resolution 
requirements (e.g. in which regions of the computational domain the smallest scale features will 
develop) are not known a priori. 

In numerical relativity, the use of AMR has played a crucial role in the study of black hole 
critical phenomena, one class of which is characterized by self-similar solutions, which thus involve 
features on arbitrarily small scales (see [7] for example). Other examples where AMR is likely to 
ultimately play an important role include the calculations aimed at predicting the gravitational 
waves generated from the inspiral and merger of a black hole binary. Here one basic length scale 
is set by the mass, M, of one of the holes (we assume roughly equal mass black holes) — clearly, 
in the vicinity of the black holes our finite difference mesh must be sufficiently fine to resolve that 
length scale. On the other hand, the gravitational waves that are generated in the late stages of 
the binary evolution will have wavelengths of order 10M, and the computational grid may have 
to extend to a distance of order 100M from the sources in order that the wave train that would 
register in a terrestrial detector can be accurately read off. Use of a single uniform grid with h 
chosen small enough to resolve the smallest scale features would not only be extremely inefficient 
in such a case, it would likely be prohibitively expensive (both in terms of CPU time and memory), 
even with the largest computers currently available. 

One can attempt to deal with the multi-scale nature of such a problem through the use of a 
non-uniform grid (or several component grids, each of which is uniform) so that the local mesh 
spacing at least roughly matches the local resolution requirement. Here we are envisioning the 
use of a priori information concerning the nature of the solution, and would further assume that 
the resolution demands are static, so that, for example, the regions requiring highest resolution 
would not move around during the course of the solution. Such a refinement strategy is often 
known as fixed mesh refinement (FMR) and has been used successfully in some recent calculations 
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in numerical relativity (see, for instance, [491 150j). However, in many situations — such as the 
modeling of objects that propagate through the computational domain — we will not have detailed 
a priori information concerning the resolution requirements, and thus algorithms that adaptively 
respond to the local development of solution features become attractive. Note that such algorithms 
need to have the capability both to introduce additional resolution when and where it is needed and 
to coarsen the local discretization scale when and where resolution demands become less stringent. 

Here wc briefly describe the version of AMR introduced by Berger and Oliger in 1984 6 , which 
has been previously extended to accommodate the treatment of elliptic equations and then incor- 
porated into the axisymmetric code, graxi [HI 03 03 E2] , used to generate the results discussed 
in Chap. 03 We note that the graxi implementation of AMR introduces several simplifying re- 
strictions relative to the original Berger and Oliger work — the most important of these include the 
requirements that (1) the computational domain be covered by a single base- level grid, (2) all com- 
ponent grids be aligned with the coordinate axes, (3) each level, I, of refinement be characterized 
by a unique spatial discretization scale, hi, (4) the refinement factor between any two successive 
discretization levels be a unique integer p, and (5) distinct grids at the same level are not allowed 
to overlap. 

The Berger and Oliger algorithm was designed for the finite difference solution of general systems 
of hyperbolic partial differential equations. Adaptivity is achieved via the use of a hierarchy of 
component grids, where each component has uniform mesh spacings in each of the coordinate 
directions (and in our case, the same mesh spacing in all directions). Specifically, the hierarchy 
is characterized by a sequence of discretization levels, I = 0, 1,2, • • ■l Ta3X ,, where I — represents 
the coarsest level of grid, and Z max the finest. Associated with each level is a spatial discretization 
scale, hi, and the scales at different levels satisfy 



where p > 2 is an integer (typically 2 in our calculations). The grid hierarchy can be denoted G;.j, 
where the index I ranges over the discretization levels, while the index j, for each level, ranges over 
the number of component grids at that level. At the coarsest level, / = 0, there is only a single 
grid, Go,0) which covers the entire computational domain. At finer levels, I > 1, there can be an 
arbitrary number of grids each of which is required to (1) have extrema (x m i n , £ max , y m in etc.) that 
coincide with mesh points of some level I — 1 grid, and (2) be completely contained within that 
grid. The level I — 1 grid, Gi-ij> that contains Gij is called the parent grid, while Gi,j is known 
as a child of Gi-\.y . As mentioned previously, distinct grids at the same level are not allowed to 
overlap. 

Fig. 13.31 shows a schematic representation of a typical grid hierarchy for a 1+1 dimensional 
problem with Z max = 2 and a refinement ratio p — 2. Note that the mesh refinement occurs in time 
as well as in space, and that the discrete time steps, ki, used in the various levels of the hierarchy 
satisfy 



hi = phi + i , 



(3.22) 



h = pki +1 , 



(3.23) 



(i.e. a constant Courant factor is maintained across the hierarchy). Also notice that the structure 
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of the hierarchy changes in time — resolution is increased where needed through the introduction of 
new grids, and decreased where it is not needed through the deletion of one or more existing grids. 
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Figure 3.3: Schematic representation of a typical grid structure for a 1+1 calculation using Berger 
and Oliger adaptive mesh refinement (AMR). In this figure we have i max = 2 and 



The core of the Berger and Oliger algorithm is a recursive time stepping algorithm whereby 
equations on coarse grids are advanced before those on fine grids. This allows the fine grid values 
at so-called refinement boundaries (i.e. boundaries of grids Gi.j,l > that do not coincide with 
boundaries of the global computational domain) to be computed via interpolation of parental values. 
Otherwise the same set of finite difference equations (including discretized versions of the physical 
boundary conditions) are used to update discrete unknowns on all component grids. Once a single 
time step has been completed at level I, and if any level I + 1 grids exist, the recursive time stepping 
algorithm is invoked to take p steps at level I + 1, so that the fine grid equations are integrated to 
the same physical time as the coarse ones. 

Another key component of the AMR algorithm is the regridding procedure in which resolution 
demands are periodically assessed, and the grid hierarchy is correspondingly reconfigured. The 
crucial task here is to decide when and where a particular solution region requires refinement or 
coarsening. Although other approaches are possible, graxi implements the strategy described in [B| 
that ties regridding to local (truncation) error estimates based on Richardson- type procedures. To 
illustrate this method, assume that our FDA can be written as 

u n+l = Qh u n ^ ( 3 _ 24 ) 

where u n+1 and u n are the advanced (t = t n+1 ) and retarded (t = t n ) finite difference unknowns, 
respectively (spatial grid indices are suppressed), and Q h is the finite difference update operator 
that advances the solution from one discrete time to the next. We can then generate a local estimate 
of the solution error by taking two time steps on the base mesh, characterized by discretization 
scale h, and then comparing the result to what we get by advancing the retarded data by one step 
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on a coarsened mesh with discretization scale 2h. That is we compute 

e h = (Q h Q h - Q 2h ) u n , (3.25) 

where the coarse grid update operator Q 2h employs the same FDA as Q h . (This technique is di- 
rectly analogous to the procedure commonly used in adaptive ordinary differential equation (ODE) 
integrators.) This error estimation can be performed dynamically and makes use of the same FDA 
that is used in the basic time stepping procedure previously described. We also note that the tech- 
nique can be applied independently of the particular PDE being solved or the specific difference 
technique that is being used. 

Once the error estimate is calculated, we can determine those grid points where the error 
estimate exceeds some predefined tolerance. These points are then organized into clusters and the 
regridding procedure is carried out. This includes creation of a new grid hierarchy (or some portion 
of it), transfer of values from the old hierarchy to the new, and interpolation of grid function values 
from parental to child grids for those spatial locations where a given level of refinement did not 
previously exist. 

A final key aspect of the Berger and Oliger algorithm involves the fine-to-coarse transfer of 
values from level I child grids to their level I — 1 parents once the level I integration has been 
advanced to the time of the level I — 1 solution. Were this not done, the level I — 1 values coinciding 
with refined regions — which by definition are not of sufficient accuracy to satisfy the error control 
criterion — would eventually "pollute" level / — 1 values in regions that should, in principle, not 
require refinement. 

Again, the reader who is interested in further details of this AMR algorithm can consult [HJ |S] . 

3.4 Excision Techniques 

Starting in the 1960s, R. Penrose and S. Hawking proved a series of theorems which are now 
collectively called the singularity theorems. The essential implication of these theorems is that 
physical singularities are generic features of spacetime, and are not a consequence, for example, of 
symmetries that a specific spacetime might possess. In particular, these mathematically rigorous 
results, as well as the results from detailed numerical simulation of gravitational collapse make 
it clear that black holes invariably contain physical singularities. This fact has proven to be 
of enormous significance for the numerical simulation of black hole spacetimes. At singularities 
physical quantities become undefined and it is meaningless to evolve spacetime there. 2 It thus 
seems to be an unavoidable conclusion that numerical simulations must somehow avoid black hole 
singularities. The traditional approach in numerical relativity (used in virtually all calculations 
prior to the early 1990's) was to use coordinate freedom to accomplish singularity avoidance. 
For instance the maximal slicing condition (described in more detail in Chap. UJ, whereby the 
trace of the extrinsic curvature tensor is required to vanish at all times, generates slices that 
automatically "slow down time" in the vicinity of a black hole singularity, essentially freezing 

2 In fact the definition of spacetime should not contain singular events. 
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the evolution before the singularity is encountered. Unfortunately, singularity-avoiding coordinate 
systems generically become pathological (i.e. coordinate pathologies develop) on dynamical time 
scales, and this typically leads to a breakdown of simulations (i.e. the code "crashes" due to an 
inability to resolve extremely steep gradients etc.). Thus the use of singularity avoiding coordinates 
alone does not appear to provide an effective resolution to the problem of long-time evolution of 
black hole spacetimes. 

An alternate approach to singularity avoidance that has become quite popular over the past 
decade or so is based on an idea due to Unruh (see for a discussion of the method). Unruh 
argued in the early 1980's that since the interiors of black holes are causally disconnected from the 
exterior universe, it is not necessary (and is in fact wasteful) to evolve regions inside black holes 
when our main interest is in the exterior region of spacetime. He thus proposed that one excise the 
interiors of black holes from the computational domain. However, since the location of a black hole 
surface — the event horizon — requires knowledge of the complete spacetime, one cannot in principle 
excise a black hole until the global solution is in hand. Unruh's second suggestion was thus to 
use the location of an apparent horizon (which one can very loosely view as an instantaneous 
approximation to the intersection of an event horizon with a given spacelike hypersurface) as the 
excision surface. Black hole excision was first explored numerically by Thornburg |53| . while the 
first successful implementation of the method in a dynamical context was due to Scidcl and Sucn 

. Some other recent references that discuss excision include [23 1523 EH EM EH1 1523 El E21 EH! • 
We will describe this method in more detail in section IP1 



3.5 Spatial Compact ificat ion 

In this section we discuss issues related to compactification of a spatial domain. This technique is 
used in the construction of the stationary boson star solutions described in Sec. 14. 3. 31 and Sec. 15.21 
For simplicity of presentation we assume here that the domain to be compactified is 1-dimensional. 

The principal motivating factor for using a compactified spatial domain is the desire to use 
exact, rather than approximate, boundary conditions in the solution of our PDEs. Specifically, 
for the problems considered in Sec. 14.3.31 and Sec. 15.21 we know the exact boundary conditions 
at spatial infinity, r = oo, whereas we generally know only the asymptotic behaviour of given 
functions at any specific finite distance, r max . Moreover, even if we know the exact form (e.g. 
precise falloff conditions) for our functions, it may be computationally advantageous to work with 
Dirichlet conditions (as the conditions at infinity tend to be), rather than the mixed (Robin) 
boundary conditions that typically result from the use of known falloff behaviour. We thus would 
like to explore one class of transformations that allows us to perform compactification. 

Suppose we want to compactify r € [0, oo) to £ £ [a, b] via the following transformation: 

where f(r) is some arbitrary function of r. To ensure C(0) = a and C(°°) — ^ we have 
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/(0) - z , (3.27) 

1 — a 

f(oo) = (3.28) 



We also wish the coordinate transformation to be regular, so that 

d( r 



dr (1 + ff 

or simply 



^0, (3.29) 



/V0. (3.30) 

In other words, / should be a strictly increasing (or decreasing) function of r. I|3.27|l . (|3.28|l and 
(|3.3Q(I are the only restrictions on a general compactifying function /. However, as we will see in 
Sec. 15. 21 when we consider the generation of initial data representing rotating boson stars, once we 
have fixed coordinates, we find that as the angular momentum of the star increases, the location 
of the extremum of the solution moves towards a boundary of the computational domain. In 
addition, for larger angular momentum values, the solution function becomes increasingly peaked, 
which means that for fixed solution accuracy and coordinates, finer resolution is required as the 
angular momentum increases. Thus, in addition to compactifying the domain, we attempt to map 
the location of the function extremum to the central region of the computational domain. 

We therefore impose the additional conditions that the coordinate transformation map some 
arbitrary points, r^, to some given values Q. We then have 

/(n) = I ^7, (3.31) 
where i = 0, • • • , N. If we assume b = 1, then one such possibility is to write 

f(r) = a + air + a 2 r 2 H 

i 

whereby can be chosen to satisfy conditions (|3.27|) and (|3.31|) . More specifically, if we want to 
compactify the radial coordinate to [0, 1], and we want to bring one particular point r = ro (say 
the location where the grid function is a maximum) to £ = 1/2, then qq = 0, aj. = 1/ro and 
ct2 = &3 = ■ ■ ■ = 0', i.e. we have 



CM 



1 + air 
r 



ro 



r 



(3.33) 
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In other words, transformation (|3.33|l maps [0, oo) to [0, 1] with C( r o) = \- 
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CHAPTER 4 

BOSON STARS IN SPHERICAL SYMMETRY 

In this chapter we study boson stars in spherical symmetry. In spherically symmetric spacetime the 
equations of motions are greatly simplified, and with a proper choice of coordinates the number of 
variables that must be evolved is significantly reduced. Therefore it is relatively easy to study the 
system numerically compared to higher dimensional (2-dimensional and 3-dimensional) simulations. 
At the same time, of all of the symmetries that could be imposed to reduce the field equations to a 
set of PDEs in one space dimension and time, spherical symmetry is clearly the most appropriate 
for the study of isolated, gravitationally compact objects. 

This chapter consists of two parts. The first part concerns the construction of static solutions 
which represent spherically symmetric boson stars in the ground state. These solutions provide the 
initial data for the simulations of critical phenomena of boson stars in spherical symmetry studied 
in the second part of this chapter, as well as for the simulation of critical phenomena of boson stars 
in axisymmetry in Chap. Together with rotating solutions computed in the Newtonian limit, 
the solutions also provide initial guesses in the construction of the general relativistic, stationary, 
rotating boson star solutions which will also be described in Chap. |S| 

The second part of the chapter concerns the dynamics of spherical boson stars. More specifically, 
in an extension of the work performed by Hawley and Choptuik [2] , we study the critical phenomena 
of boson stars which are driven to the threshold of black hole formation via an external perturbing 
agent. As in the previous work, our boson stars are perturbed by a real massless scalar field 
and an overall amplitude of this real scalar field is used to tune to criticality. However, we use 
a different coordinate system — maximal- isotropic coordinates — than that used in 0. Although 
the accuracy of finite difference calculations in any given coordinate system can in principle be 
estimated using intrinsic means (e.g. convergence tests), we feel that it is nonetheless useful to 
reproduce the calculations of P3 G] in a different coordinate system. In addition, since maximal- 
isotropic coordinates can penetrate black hole horizons (the polar-areal coordinates used in [2] 
generally cannot), we can also implement black hole excision in our simulations. We emphasize 
however, that we feel that our excision calculations — which are restricted to spherical symmetry 
and use an ad hoc boundary condition at the excision surface for the lapse — are only a minor result 
of the thesis. 

The principal new results presented in this chapter concern the long-time evolutions of subcriti- 
cal simulations in critical collapse of boson stars. We provide numerical evidence that — contrary to 
the previous conjecture that subcritical solutions disperse most of the original mass of the boson 
star to large distances — the late time behaviour of subcritical evolution is characterized by oscilla- 
tion about a stable boson star solution. We also apply a linear perturbation analysis similar to that 
in j2| and confirm that the observed oscillation modes agree with the fundamental modes given by 
perturbation theory. (We use a code provided by S. Hawley to generate the frequencies from 
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the perturbation analysis.) We also observe that there is generically an overall lower-frequency 
modulation of the post-critical-phase oscillations. At this time we have no explanation for the 
origin of this additional oscillatory behaviour, although we can safely rule out the possibility that 
it originates from the beating between the fundamental and first harmonic mode. 

The outline of the remainder of the chapter is as follows. In Sec. 14.11 we derive the equations 
of motion for boson stars in a general spherically symmetric spacetime. In Sec. 14.21 we specialize 
the equations of motion to maximal-isotropic coordinates. In Sec. 14.31 we derive systems of equa- 
tions for the initial value problem in maximal-isotropic coordinates, polar-areal coordinates and 
compactified maximal-isotropic coordinates, and briefly describe the family of stationary solutions 
that we can generate using any of the three approaches. In Sec. 14.41 we present the results of near- 
critical evolutions and, for the case of subcritical evolutions, compare the oscillatory behaviour 
observed at late times to that predicted from perturbation theory. Finally in Sec. 14. 51 we study the 
implementation of black hole excision technique applied to our system. Throughout this chapter 
we denote ' = d r and '= d t , and m = 1 is used for the particle mass of the boson field. 



4.1 Spherically Symmetric Spacetime 

Loosely speaking, we may regard a spherically symmetric spacetime as one that admits a preferred 
timelike observer such that the spacetime is spherically symmetric about every point on this spe- 
cial observer's world-line. 1 The most general metric for a time-dependent spherically-symmetric 
spacetime can be written as jHSI H20 EH] 

ds 2 = (-a 2 + a 2 (3 2 ) dt 2 + 2a 2 f3dtdr + a 2 dr 2 + r 2 b 2 dtf , (4.1) 

where dfl 2 = d9 2 + sin 2 9 dtp 2 is the metric on the unit 2-sphere, (3 is the r-component of the shift 
vector, (3 l = (/3, 0, 0), and a, /?, a, b are functions of t and r only. Corresponding to the above metric 
is an extrinsic curvature tensor with only two independent components 2 

K*j = dizg (K r r ,K%,K%) , (4.2) 

where K r r , K e g are also functions of t and r. Together with a, /?, a and b we thus must deal with a 
maximum of 6 out of 16 possible geometrical variables in spherically symmetric calculations. The 
4-metric can be written in matrix form as 

" -a 2 + a 2 p 2 a 2 (3 

a 2 (3 a 2 

■ hl "~ (rb) 2 

(r&sin(9) 2 

1 Spherically symmetric spacetime is defined as one which admits the group SO (3) as a group of isometries, with 
the group orbits spacelike two-surfaces. 

2 The form of the extrinsic curvature can be inferred by considering the expression dfjij = — 2a / yn t K k j +/3 k ^ij t k + 



(4.3) 
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from which it immediately follows that the inverse metric is given by 



-1/a 2 (3/a 2 

[i/a 2 l/a 2 ~/3 2 /a 2 

l/(rb) 2 

l/(r& sin <9) 2 



(4.4) 



We also have \/ = g = aab 2 r 2 sin 0, where g is the determinant of the 4-metric. Now considering the 
3- metric, 7^, we have jij — gij, so 



a 2 
(rb) 2 
(rb sinO) 2 



and 



7*-? = 



1/a 2 
l/(rb) 2 
l/(r&sin6>) 2 _ 

The non-vanishing Christoffel symbols constructed from 7^ are 



1 r r 


a' 
a 


T r 99 


(r 2 6 2 )' 
2a 2 ' 


rr 4><t> 


= sin 2 er r 99, 


r e 

1 r9 


r e (rb)' 

= J- Or = , 

ro 


r e 

l 4,4, 


= — sin cos 9 , 




= r = T 8 r e 




= r% e = C0tfl 



(4.5) 



(4.6) 



(4.7) 



The non- vanishing components of the Ricci tensor R l j are 

■(rb)' ' 



R T r — 



R% = 



2 
arb 

1 

R 6 e. 



rb 



rb)' 



(4.8) 

(4.9) 
(4.10) 
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and the scalar curvature R is 
R = 



R r r + R t 

R r r + 2R e 



2 


~(rb)'- 




2 


arb 




a 




a(rb) 2 


2 




\(rb)>- 


' 1 


arb 




a 




+ -T 

rb 



rb 



(rb)' 



-(rb)') -a 
a 



(4.11) 



On the other hand, from (|2.44() . (|2.11() l|2.12[) . (|2.18|) and l|2.19|l the non- vanishing components of 
the matter source terms are 



P 



Jr = 



s r , 



m 2 + m 2 __ 

2a 2 2 ' 

n*$ + n$* 2 ._ 

2a— = a f ' 

p-m 2 ^ 2 , 



2a 2 



= S" 



S 



3|n| 2 H$| 

2a 2 



-m 2 \(j)\' 



Here, we have defined the auxiliary scalar-field variables 



$ = </>' 
a 
a 



n 



(4.12) 

(4.13) 
(4.14) 
(4.15) 
(4.16) 
(4.17) 



(4.18) 
(4.19) 



which are also useful in recasting the Klein-Gordon equation in first-order-in-time form. 
The Hamiltonian constraint (|2.13|) now becomes 





IW" 


' 1 


arb 1 


a 


+ -T 

rb 



rb 



rb)' 



^K r r K% + 2K 



o 



8tt 



l$l 2 + ini 



(4.20) 



while the momentum constraint (|2.14|) is 

, (rb)' 



K° 



rb 



(K 6 



2-7T 

-k\) = — (ir$ + n$*) 

' a 



(4.21) 
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A straightforward but tedious calculation of the evolution equations gives 
a = —aaK r r + (a/3)' , 



b = -abK g + ^-(rby 



K r r 



K° 



l3K r ' r 



1 (ct 



a 





Urb)n 


arb 


a 



KK r 



An 



(rb) 2 a(rb) 2 



arb / 

( rb ) 

a 



a(KK e e -ATrm 2 \(j)\ 2 ^ 



2|#| 



(4.22) 
(4.23) 

(4.24) 
(4.25) 



The definitions (|4 . 1 8J) . Ij4.19|) and the Klein-Gordon equation (|2.45|1 leads to the evolution equations 
for the scalar field and the auxiliary field 



-n + /3$. 



ri 



(rb) 2 (/3IL + 
V a 



(rb) 



aam cj) + 2 



(r&y 
rb 



n. 



(4.26) 
(4.27) 
(4.28) 



Equations H4.2Ufl - H4.28fl comprise the general system of equations for a general-relativistic complex 
scalar field in spherical symmetry. To study the system numerically, one needs to choose specific 
coordinate conditions, impose boundary conditions and set up appropriate initial conditions. In the 
following section we will specialize the above equations to the case of maximal-isotropic coordinates. 



4.2 Maximal-isotropic Coordinates 
4.2.1 The equations of motion 

The equations presented in the previous section are valid in any spherically symmetric coordinate 
system compatible with (|4.1f) . However, in order to perform numerical simulations we need to 
choose a specific coordinate system, which, in the 3+1 approach is equivalent to prescribing the 
lapse, a, and the shift vector component, f3. For a variety of reasons, we have chosen to adopt 
so-called maximal-isotropic coordinates. First, these coordinates have been used in several pre- 
vious calculations in numerical relativity (for example, see |(i8l and have generally worked 
well. Second, maximal-isotropic coordinates are the specialization to spherical symmetry of the 
coordinate system used in the axisymmetric code discussed in the next chapter. This fact allows 
us, in principle, to directly compare results from the spherical code to those from the axisymmetric 
code, when the latter is supplied with spherically symmetric initial data. Third, in contrast to 
the polar- areal coordinates used in maximal-isotropic coordinates can penetrate apparent 

horizons, and thus using them, we can more readily study the process of black hole formation. In 
particular we can incorporate black hole excising techniques as discussed in Sec. 14.51 

The "maximal" part of maximal-isotropic refers to the maximal slicing condition, which fixes 
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the lapse function. A hypersurface is maximal if its mean extrinsic curvature vanishes; that is if 

K = K\ = . (4.29) 

Geometrically, it can be shown that K — implies that the volume of the hypersurface is maximized 
with respect to localized, infinitesimal deformations of the slice (analogous to the minimization of 
the surface area by a soap film bounded by a wire frame). Computationally, we implement maximal 
slicing by choosing initial data for the extrinsic curvature that satisfies K = 0, and then demanding 
that 

K(t,r) = 0. (4.30) 

for all t and r. We note that the maximal condition allows us to eliminate one of the two non-trivial 
components of the extrinsic curvature from the computational scheme. In particular, since 



K = 



1 



-K r 



(4.31) 



we choose to eliminate K e g. As we will see shortly, the maximal condition leads to a second order 
ODE that constrains the lapse function, a at all times. 

Given maximal slicing, the remaining (spatial) coordinate freedom is fixed by the isotropic 
condition which requires 

a = b=iP(t,r) 2 , (4.32) 

for some positive function ip(t, r), such that the spatial metric takes the "isotropic" (or 3-conformally 
flat form) 

^ds 2 =^(dr 2 +r 2 <m 2 ) . (4.33) 
This is implemented by choosing initial data such that a(0, r) = &(0, r) and then demanding that 

a(t,r) = b(t,r), (4.34) 

for all t and r. As will be shown below, this condition leads to a first-order ODE that fixes the 
shift vector component, /3, on each slice. 

Turning now to the equation for the lapse, a, that is implied by K = 0, we note that by taking 
the 3-trace of (I2.2UII we have: 



K = (IK' - 



1 



a{rbf 



{rbf ; 
a 



a 



[R + K 2 + 4tt (S - 3p)] 



(4.35) 



Using K = 0, as well as l|2.13|l . (|4.17|l and l|4.12|l . for R, S and p respectively, we have 



or) 




a' + 







Ana 2 m 2 \<t>\ 2 - 87r|n| 2 - -a 2 K\ 2 



a = 0. 



(4.36) 
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Using the isotropic condition, a = b = ip 2 , this can be further simplified to 
[(V-r) 2 ]' , 



a" + l ( 1 < 

(jpr) 2 



47r^ 4 m 2 |0| 2 - 87r|n| 2 - | (ifj 2 K r r 



a = 0, (4.37) 



rip 2 dr 2 ^ ^ 



An^m 2 \4>\ 2 - 8^|II| 2 - | {^ 2 K r r )' 



a = , (4.38) 



where d/dr 2 is the derivative with respect to r 2 . We note that here and elsewhere, and following 
Evans [7Dj, we often rewrite terms of the form d r f(t,r) as pr p ~ 1 d r p f(r,t) for functions / satisfy- 
ing lim r _,o f{r, t) — r p f p (t) + 0(r p+2 ). This technique ensures that when standard second-order 
centered difference formulae are applied to such terms, the leading order regularity behaviour is 
preserved in the discrete domain as r — > 0. 

The ODE that fixes the shift vector component j3 (i.e. the "isotropic condition for the shift"), 
can be easily derived from (|2.17|) . by equating the right hand sides of the respective evolution 
equations for a(t,r) and bit, r). Doing this we find 

r(?\ =a {K% ~ K%) , (4.39) 
which can be further simplified using the maximal condition to yield 

r (r) = \ aKr -- ( 4 - 4 °) 

Both the slicing condition (|4.38|l and the shift-component equation (|4.4U|) must be supplemented 
by appropriate boundary conditions, but we will relegate these details to App. 

Having fixed the kinematical geometric variables (i.e. a and /?), we now consider the specific 
form of the equations of motion for the dynamical geometric variables, as well as those for the 
scalar fields that are coupled to the gravitational field. 

Due to our restriction to spherical symmetry (in which case the gravitational field has no 
dynamics that is not tied to the dynamics of a matter field), as well as our choice of maximal- 
isotropic coordinates we can implement a so-called fully constrained scheme that uses only the 
constraint equations to update the dynamical geometric variables in time. Given our coordinate 
choice, the only non-trivial geometric variables that remain are the (3-)conformal factor, ip, and 
the extrinsic curvature component, K r r . These can be determined from the Hamiltonian and 
momentum constraints, respectively. 

Specifically, using expressions lj!H3|) . gZD} , CED), (O^ - 6321), lETT^ and <|4?T5)l . we have 

3 d ( ^\ . K, r , _ ^l*l 2 + l n l! + TOW V (4 . 41) 



+ — K r 2 - 
+ 16 r 


-< 


12 K r = 


An 



ip 5 dr 3 \ dr J 16 \ V' 4 

(rilj 2 )' 47r 
K r r ' + r r = -— (n*$ + n$*). (4.42) 



(|4.38|) . 14.4011 . (|4.41|) and (|4.42|) completely determine the geometric variables a,j3,ip and K r r in 
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maximal-isotropic coordinates. 

A caveat is in order here. A fully-constrained evolution works well as long as we are not 
implementing black hole excision. When excision is being performed, we must supply boundary 
conditions for ip and K r r at the excision surface, and these are provided by the evolution equations 
for those variables, which are 



• i (^py 
i> = - 5 o^r r + — 



(3K T 



2a 



[rip 



r 2 ip e 



inm 2 a\(j)\ z 



(4.43) 
(4.44) 



Now considering the complex-scalar matter field, and using l|4.26JI . (|4.27|) . 14.28fl and i|4.32[l . we 
have the following equations of motion: 



a 

?2 



$ = 



n 



3 d 
^dr 1 



r^ 4 I /3n + — $ 



aip m 4> - 



(rip 2 )' 
rip 2 



aK\ + 2/3 



n. 



(4.45) 
(4.46) 
(4.47) 



()4.38jl . H4.40jl . (|4.41jl . (|4.42|) . I|4.45j) - I|4.47j) constitute the basic set of equations that we use in our 
study of boson star dynamics in maximal-isotropic coordinates. 

As discussed previously, in order to gravitationally perturb boson stars, we will incorporate 
an additional, minimally coupled, massless real scalar field in our model. It is straightforward to 
generalize the above equations of motion to account for such an addition, and here we will merely 
state the results. To simplify notation we denote the real and imaginary parts of the complex scalar 
field as (pi(t,r) and (p2(t,r), respectively, while (p3(t,r) is the massless real scalar field: 



rip 2 dr 2 ^ 



4n^m 2 ]T <P 2 - Stt - | {^K\f 



a = 0. 



(4.48) 



(4.49) 
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(4.50) 
(4.51) 



i=l 



CHAPTER 4. BOSON STARS IN SPHERICAL SYMMETRY 



42 



n, 



-2^ + /»< 



(4.52) 
(4.53) 



3 d 

aK r 



"1 0i (1 — <5 i3 ) 



2/3 



(4.54) 



4.2.2 The mass aspect function 



Before discussing the construction of initial data, we want to note that for diagnostic purposes it 
is often useful to compute the mass aspect function M(t,r), which approaches the ADM mass as 
r — ► oo. In Schwarzschild-likc coordinates (T, R, 9, <fi), (b = 1, (3 = in l|4.1|) ) the mass aspect can 
be easily expressed in terms of the metric function a, since for large R (i.e. outside the support of 
any matter fields in the spacetime) we expect the metric to be of the form 



ds 2 = - 1 - 



2M 
~R 



dT 2 



1 - 



2M 
~R 



dR 2 + R 2 dVL 2 



(4.55) 



where M is a constant to be interpreted as the (ADM) mass of the system. On the other hand, 
the metric (|4.1|) specialized to the case b = 1, (3 = is 



ds 2 = -u 2 dT 2 + a 2 dR 2 + R 2 dil 2 . 
We can thus define the mass aspect function M(T, R) such that 

2M (T, R) N 



1 - 



R 



or, solving for M, 



M(T,R)=%[1 



= a(T, R) 1 



1 

a(T,Rf 



(4.56) 



(4.57) 



(4.58) 



We can generalize this expression to the case of an arbitrary time-dependent spherically-symmetric 
coordinate system as follows. We first note that, in spherical symmetry, and at least in vacuum 
regions, M is a geometric (i.e. coordinate-independent) quantity, as is the areal radius R (note that 
we generally have R(t,r) = rb(t,r) from l|I3J). 

We now consider the square of the gradient of the areal radius, V A1 i?V A 'i?, which is also a 
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geometrical invariant. In Schwarzschild coordinates we have 

= g RR V R RV R R 

= g RR 

2M 
= '-If 

Thus we have the coordinate-independent formula 



(4.59) 



M = — (1 — V p i?Vfl) , (4.60) 



where, again, we emphasize that R is the areal radius. 
In maximal-isotropic coordinates we have R — ip 2 r, so 



M[t, r) = tJL [i - V M (V 2 r) V (^ 2 r)] . (4.61) 
With a little algebra the above expression can be written as 

M(t, r) = {^f^j K\ 2 ~ 2A 2 (V> + rtl)') . (4.62) 

Finally, we note that we summarize all essential equations, boundary conditions, as well as the 
finite differencing of the equations discussed in this section in App. [5] 

4.3 The Initial Value Problem in Spherical Symmetry 

As we have discussed previously, we are particularly interested in the study of the dynamics of boson 
stars which are stationary (or equivalently static in the current case of spherical symmetry) solutions 
of the coupled Einstein-Klein-Gordon equations. Determination of initial data for the scalar and 
gravitational fields that represents such stars is a special case of the general initial value problem 
for our model. In fact, we have used three different methods to obtain spherically-symmetric 
boson star initial data. The three approaches are summarized in the following sub-sections: each 
is conveniently labelled by the spatial coordinates used in the initial value computation, and each 
has its own specific advantages. 



4.3.1 Maximal-isotropic Coordinates 

We first discuss the construction of boson star data within the context of the maximal-isotropic 
coordinate system described in the previous section. To this point we have been considering the full 
dynamical equations describing the Klein-Gordon field coupled to gravity in spherical symmetry. 
To construct a star-like configuration we need to impose additional specific constraints on the scalar 
field. Ideally one would like a "star" to be described by a localized, time-independent matter source 
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that generates an everywhere regular (i.e. non-singular) gravitational field. However, for the case 
of a complex scalar field, it can be shown that such regular, time-independent configurations do 
not exist [7T]. Despite this fact, since the stress-energy tensor ()2.44|) depends only on the modulus 
of the scalar field (and the gradients of the modulus), one can construct scalar field configurations 
with harmonic time-dependence that produce time-independent metrics. Specifically, we adopt the 
following ansatz for boson stars in spherical symmetry: 

<j>(t,r)=Mr)e- iujt , (4.63) 

and then demand that the spacetime be static, i.e. we demand that the metric admits a timelike 
Killing vector field \ which is orthogonal to the t — const, surfaces. Adapting coordinate time to 
the timelike Killing vector field, we have 

(3 = 0, (4.64) 



for all time t. Additionally, we have that the time derivatives of any of the geometrical variables 
identically vanish. It then follows immediately from ((4.43(1 that 

K r r = 0. (4.65) 

As is necessary for the consistency of the ansatz (I4.63|) . the isotropic condition for (3 (14.40(1 is 
automatically satisfied, and we are left with geometrical variables a(0, r), ip(0,r) and <fio(r) that 
need to be determined from the maximal slicing condition 1(4.38(1 , the Hamiltonian constraint ((4.41(1 
and the Klein-Gordon equation 1(4.47(1 . respectively. Before considering the solution of those three 
equations, we note that from the ansatz 1(4.63(1 we have 

<j> = -iufo e~ luJt , (4.66) 

(/>' = 0o'e- 4 "*, (4.67) 

and hence 

ib 2 

II = -i— oj(j) e- %ut , (4.68) 
a 

$ = o 'e- iwt . (4.69) 

Therefore, 

lT$ + n$*=0, (4.70) 

and the momentum constraint ((4.42|) is satisfied for all times. 

Using the ansatz ((4.631) for t = 0, as well as (3 = K r r = 0, straightforward manipulation of the 
Hamiltonian constraint, the Klein-Gordon equation, and the maximal slicing equation yields the 
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following set of ODEs: 



2* 



2 I J.2 



= A , 



— I h — $ + -0 m 

r a y / 



r yj J \ cr 



(4.71) 
(4.72) 
(4.73) 
(4.74) 
(4.75) 
(4.76) 



Here, in order to simplify notation, we have dropped the subscript "0" , making the identifications 
</>(r) = 4>o( r ) and $(r) = 4>'(r) = <j>' (r). We have also introduced auxiliary variables ^f(r) = ip'{ r )> 
$(r) = 4>'{r) and A{r) = a'(r) in order to cast the above system of nonlinear ODEs in a canonical 
first-order form. We assert that for any given value of 0(0) = 4>o(0), the system (|4.71f) - (|4.76f) 
constitutes an eigenvalue problem with eigenvalue lo = w(</>(0)). That is, for any specific value of 
(j)(Q) (which one can loosely view as being related to the central density of the star), a solution 
of (|4.63|l that satisfies the appropriate regularity and boundary conditions will only exist for some 
specific value of ui. 3 The system H4.7Hl - l|4.76|l must be supplemented by boundary conditions, 
some of which are naturally applied at r = 0, with the rest naturally set at r = oo. In particular, 
regularity at r — implies 



*(0) 
$(0) 

A(0) 



0. 

o, 

0. 



while at the outer boundary, we have 



lim ip(r) = 1 — 



lim Mr) 

r — 'too 

lim a(r) 



0, 
2 



C 



1 



(4.77) 
(4.78) 
(4.79) 



(4.80) 
(4.81) 

(4.82) 
0. 



Here the second condition follows from the expectation that (j) should decay exponentially as r 
(The outer boundary conditions for ip and a are derived in detail in Add, loll 

We further note that due to the homogeneity and linearity of the slicing equation, we can always 

3 In fact, in general one expects an infinite, discrete set of eigenvalues, u>i, i = 0, 1,, • • • (for any value of </>(0)), 
corresponding to "wave functions" 4>i(r) with increasing number of nodes. We restrict attention here to "ground 
state" boson stars, where <j>o{r) has no nodes. 
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arbitrarily (and conveniently) choose the central value of the lapse via 

a(0) = 1 , (4.83) 

and then, after integration of (|4.71|) - (|4.7tj|) . can rescale a and uj simultaneously to satisfy the outer 
boundary condition for a: 

a{r) — ► ca{r), (4.84) 
uj — ► cuj. (4.85) 

where c is given by 

c = g/%4= 1 , (4.86) 

and r max is the radial coordinate of the outer boundary of the computational domain. 

As mentioned above, any solution of (|4.T1|) - (|4.76|1 can be conveniently labelled by the central 
value of the modulus of the scalar field, <po(0) = 0(0). For any given value of (j>o(0), we must 
then determine the eigenvalue, u>, and in the current case of maximal-isotropic coordinate, the 
central value of the conformal factor ij>(0), so that all of the boundary conditions are satisfied. In 
principle, we can compute pairs [^,"0(0)] as a function of 0o(O) using a two-parameter "shooting" 
technique |4"5ll?2] . 

The process outlined above is the most direct approach for finding boson star solutions in 
maximal-isotropic coordinates. The principal disadvantage lies in the fact that we found it difficult 
to implement a convergent 2-parameter shooting method that allows one to compute values [uj, if>(0)] 
given O (0). 



4.3.2 Polar-Areal Coordinates 

In this subsection we describe a technique for generating boson star initial data in maximal-isotropic 
coordinates by first constructing the stars in so-called polar-areal coordinates, and then performing 
a coordinate transformation. 

Polar-areal coordinates, which have seen widespread use in spherically symmetric computations 
in numerical relativity, can be viewed as the generalization of the usual Schwarzschild coordinates 
to time- dependent, spherically symmetric spacetimes. As with maximal slicing, the slicing condition 
in this case — known as polar slicing — is expressed as a condition on the mean extrinsic curvature: 

K = K\ . (4.87) 

Since in general we have K — K l \ = K r r + 2K s g, this condition is implemented by requiring 

K e e (t,r) = K 6 g{t,r) =0, (4.88) 

for all t and r. 

The spatial coordinates are fixed by demanding that the coordinate r measure proper surface 
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area (i.e. that it be an areal coordinate), which, in terms of the general form (|4.1() implies that 

6 = 1. (4.89) 

It can be shown that this combination of polar slicing and areal spatial coordinate induces further 
simplification in the 3+1 form of the metric — namely the shift vector component (3 identically 
vanishes, so we are left with a metric of the form 



ds 2 



-a 2 dt 2 + a 2 dr 2 + r 2 dVL 2 . 



(4.90) 



As before, to construct star-like solutions, we adopt the time-harmonic ansatz (|4.63|l for the 
complex scalar field, adapt the the time coordinate to the timelike Killing vector field, and require 
the spacetime to be static. We again find that the extrinsic curvature tensor vanishes identi- 
cally (so that, for static data, the slicing is maximal as well as polar), and that the momentum 
constraint (|4.42|) is automatically satisfied. 

Again, considering the Hamiltonian constraint, the Klein-Gordon equation, and the slicing 
condition 

K% = , 



at t = 0, we have (dropping the subscript 0's as before): 
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$' = - (1 + a 2 - 47rrVmV) --(- 

' r \< 

In this case, the regularity conditions are 

o(0) = 1, 
$(0) = 0, 

while the outer boundary conditions are 

lim <f>(r) ~ , 

r — >oc 

1 

lim a(r) 



m 2 ) 6a 2 . 



a{r) 



(4.91) 

(4.92) 

(4.93) 
(4.94) 
(4.95) 



(4.96) 
(4.97) 



(4.98) 
(4.99) 



As before, we can convert the last condition to an inner condition on a by taking advantage of the 
linearity and homogeneity of the slicing equation. Specifically, we can again choose a(0) = 1, and 
then after integration of H4.92|I - I|4.95JI simultaneously rescale a(r) as well as the eigenvalue, cj, so 
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that the outer boundary condition for a is satisfied. 

We again consider the family of boson star solutions parametrized by the central value of the 
modulus of the scalar field, <f>o(0). In this case, given a value of <fro(0), and using the conditions 
a(0) = 1, a(0) = 1, $(0) = 0, we need only adjust the eigenvalue ui itself in order to generate 
a solution with the appropriate asymptotic behaviour (i.e. so that linv^oo <p(r) — 0). This is a 
classic 1-parameter shooting problem, which is comparatively easier than the 2-parameter shooting 
method described in the previous section. 

Once we have computed a solution in areal coordinates, we can perform a coordinate transfor- 
mation from areal coordinates to isotropic coordinates |73l 174] (recall that the maximal and polar 
slices coincide for the static case) . Essentially this amounts to solving an ODE of the form 

r ln=iW = 

dr 
dR 

and the details of the transformation are 
4.3.3 Compactified Maximal-Isotropic Coordinates 

The two methods discussed thus far for determining boson star initial data both use a "shooting" 
method to fix one or more parameters (an eigenvalue, and possibly an initial condition) so that the 
solution of a system of ODEs satisfies all appropriate regularity and asymptotic conditions corre- 
sponding to a star-like solution. The shooting technique is inherently iterative and, particularly 
for the two-parameter case, one needs a fairly good initial estimate of the unknown parameters for 
the iteration to converge quickly. 

An alternate method, which seems to be the best for computing boson star initial data in 
maximal-isotropic coordinates (without resort to coordinate transformation) involves the use of a 
compactified spatial coordinate and the solution of the eigenvalue problem using standard linear 
algebra software. This approach has the further advantage of generalizing quite easily to problems 
with more than one spatial dimensions; i.e. to boundary value PDEs with eigenvalues, although 
it must be stressed that naive application of standard linear algebra software for the eigenvalue 
problem is in general inefficient for higher dimensional problems (see the discussions in 15.2.20 . 
(Such inefficiency can be remedied, however, by the use of a more efficient algorithm such as the 
multigrid eigenvalue method discussed in Sec. 15. 2.151 ) Shooting methods, on the other hand, are 
generally restricted to systems of ODEs. 

The key idea here is simply to introduce a compactified coordinate, (, as previously discussed 
in Chap. |31 



I + y/ay R 
a 



J R=R n 



R 



(4.100) 



given in App. iDl 



(4.101) 
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where £ g [0, 1]. Under the above transformation we have: 

c 

V = 

i-C 

(i-C) 4 vf, 



V 2 



where V 2 ee 3^ (r 2 f), and V 2 ee 3^ (c 2 ^). 

With this definition the system of ODEs H4.71jl - I|4.7fi[l can be rewritten as 



V 2 , 
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V> 4 
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U = o. 
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H^-nrM^ = 0. 



(4.102) 

(4.103) 
(4.104) 



(4.105) 
(4.106) 
(4.107) 



As explained in Chap. the advantage of using a compactified spatial coordinate is that we can 
precisely impose the boundary conditions at infinity in our numerical computations. Specifically, 



in the current case we have </>o(C)l 



o, V(C)I 



C=i 



1 and a(C)| 



After finite differencing (|4.105|l - l|4.107|l we are left with a standard nonlinear eigenvalue problem 
that can be treated numerically using standard linear algebra packages. Details of our implemen- 
tation are described in App. El 



4.3.4 Family of Stationary Solutions 

For completeness, we display a few typical stationary solutions representing spherically symmetric 
boson stars in Fig. 14.11 Specific properties of the entire one-parameter family of stars are shown in 
Fig. 14.21 ( ADM mass vs central scalar field value </>o(0)) and Fig. 14.31 feigenvalue ui vs <j) (0)). These 
diagrams will subsequently be compared with similar plots generated for families of axisymmetric, 
rotating, stationary solutions in Chap. |SJ 

We remark that the boson stars become more compact (in other words, have smaller and smaller 
effective radii) as (f>o(0) increases. We also note that, as is apparent from Fig. 14.21 there exists a 
maximum mass, Af max 0.633Af 2 j/m, for the family of stars, corresponding to a central scalar 
field value ^o(0) ~ 0.08. Thus this model exhibits the analog of a Chandrasekhar mass limit, above 
which no static configurations exist, although we should emphasize, that in contrast to the white- 
dwarf case, the existence of an upper bound on the mass of boson stars is purely a relativistic effect. 
In particular, no such upper bound exists for non-relativistic Newtonian boson stars. Furthermore, 
perturbation analysis, as well as full simulations, show that configurations depicted in Fig. l4.2l that 
lie to the left of the mass maximum (i.e. for </>o(0) < 0.08) are stable to radial perturbations, while 
those to the right (4>o(0) > 0.08) are unstable to such perturbations. Indeed, this instability will 
play a crucial role in our analysis of critical behaviour in the model which we consider in the next 
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Figure 4.1: Static spherically symmetric boson stars. Plots of 4>o(r), tp(r), a(r) and M(r) for bo- 
son stars with </>o(0) = 0.01,0.02,0.04 and 0.07 in maximal-isotropic coordinates. In 
general, for boson stars such as these which are on the stable branch, the configu- 
rations become more compact and of higher mean density as the central scalar field 
value increases. 
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0.2 0.4 0.6 0.8 



0o(O) 

Figure 4.2: ADM mass Max>m of spherically symmetric boson stars vs the central scalar field 
value 0o (0). Note that there exists a maximum value of the ADM mass — M max w 
0.633Afpj/m — above which there are no static solutions. Note also that each ex- 
tremum point satisfy dMADM(0o(O))/rf<fo(O) = 0, which corresponds to change in 
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Figure 4.3: Eigenvalue, u>, of spherically symmetric boson stars vs the central scalar field value 
4>o (0), for the system 14.92l) - (|4.95l) . The figure on the left shows the value of log 10 (co>) 
before rescaling of a. The figure on the right shows the value of lo after rescaling. 
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section. 

4.4 Critical Phenomena of Boson Stars 

In the remainder of this chapter we focus on black-hole threshold behaviour (or critical behaviour) in 
the Einstcin-Klcin-Gordon system. In particular, following Hawley and Hawley and Choptuik [2] 
we study black-hole threshold solutions that are generated by gravitationally "perturbing" stable 
boson stars using an additional massless scalar field as the perturbing agent. The massless field 
is not explicitly coupled to the complex, boson star field, but the two fields can and do interact 
through the gravitational field. The previous calculations have provided strong evidence that 
the black hole transition in this case is Type I (so that the black hole mass just above threshold 
is generically finite) , with the critical solution itself being a perturbed boson star on the unstable 
branch (i.e. lying to the right of the mass maximum in Fig. 14. 2[) . In addition, scaling laws of the 
form 

r(p)~-7ln|p-p*|, (4.108) 

were observed for near-critical evolutions, where p is the family parameter, p* is the critical param- 
eter value that demarks the black hole threshold, r is the length of time that the evolution remains 
"close" to the (static) critical solution (i.e. the "lifetime" of the near-critical solution), and 7 is a 
scaling exponent that depends only on which of the infinitely many critical solutions is generated 
by the p = p* evolution. Furthermore, in accord with the now-familiar picture of one-mode insta- 
bility of critical solutions, the values of 7 determined from the scaling relation (|4.1U8|I were shown 
to be in agreement with the reciprocal Lyapunov exponents associated with the single unstable 
mode of the (unstable) boson star which was a best match to the mean motion of the near critical 
evolutions. In this regard we should point out that with reference to Fig. 14.21 we expect a change 
of stability at each extremum of the plot of A/adm vs (f>o(0). In particular, those configurations to 
the right of the absolute maximum of the plot, but to the left of the subsequent local minimum 
(near </>o(0) = 0.3) each have a single unstable mode in perturbation theory, and thus each can act 
as an "intermediate attractor" in critical collapse. 

As we will see, our current calculations confirm these basic results (and do so in a different 
coordinate system than that used in the previous investigation), but also provide clear indications 
that the end-state of subcritical evolutions is quite different than that which was proposed in 

4.4.1 Setup of Numerical Experiments 

The PDEs solved in the simulations discussed here are those listed in Sec. 14. 21 including the equation 
of motion for the massless scalar field. We also provide a summary of the equations of motion of 
the system, the boundary conditions, and details of the finite difference approximation used in 
Adds. IaI and IbI Additionally, results of convergence tests of the code are discussed in Add. IU1 

In order to study critical behaviour in the model we start with initial data for the complex 
field that represents a boson star on the stable branch (i.e. a star with a central scalar field value 
cpo (0) < 0.08). We generally choose a configuration that is reasonably relativistic, i.e. with <f>o(0) 
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bounded away from 0, but not too close to the instability point, (f>(0) w 0.08. 

To drive the boson star to criticality, we implode a (spherical) shell of massless scalar field on 
to it. Specifically, we choose initial data for the massless field of the following "gaussian" form 



j(0,r) = A 3 exp 



r - r 



(4.109) 



where A3, ro and a are adjustable parameters, controlling the overall amplitude, position and width, 
respectively, of the imploding gaussian wave packet. To ensure that the massless field is almost 
purely in-going at the initial time, we specify the "conjugate" variable II3 = ip 2 /a ^3 — j3<f>zj as 
follows: 

n 3 (0, r) = - ( $3(0, r) + . (4.110) 



In all of our studies described here, we have fixed ro and a in (|4.109|l to ro = 40 and a = 5. This 
ensures that the support of the massless field is well separated from that of the complex field (i.e. 
from the boson star per se) at the initial time. 

A typical evolution of initial data of the form described above proceeds as follows. Once we have 
fixed the boson star configuration, we complete the specification of the massless scalar field initial 
data by fixing the overall amplitude factor, A3, and then start the simulation. Initially, the shell 
of massless scalar field implodes towards r — at the speed of light, while the boson star "sits" 
in its static state centered at the origin. As the in-going massless shell reaches the region of space 
occupied by the boson star, its contribution to the overall gravitational field tends to compress the 
boson star to a higher mean density and smaller radius. The massless field passes through the origin 
and then "explodes" outward, eventually propagating off the computational domain. Depending 
on the strength of the perturbation from the massless field, we find that the compressed boson 
star either relaxes to something resembling a stable boson star with large-amplitude oscillations, 
or collapses to form a black hole. Thus by adjusting the massless scalar amplitude factor, A3 — 
which we generically use as the adjustable parameter, p, in our study of critical behaviour in the 
model — we can tune the evolution to the threshold of black hole formation. In practice we use a 
bisection search to refine our estimate of the critical value, and can carry the search to machine 
precision (8-byte real arithmetic), so that AA3/A3 ~ 10~ 15 . Typically we have chosen the number 
of mesh points N r — 1025 on < r < 50, a Courant factor At/Ar = 0.3, and the coefficient of 
Kreiss-Oliger dissipation Cd — 0.5 (see App. iBlfor the definition of e<j). Note that there is a whole 
family of critical solutions (see Fig. 14. 6|) for different initial data. Therefore if we, for instance, 
changed ro from ro = 40 to ro = 35 in (|4. 109(1 . we would find that A\ changes, i.e. A% — Ag(ro, a), 
and therefore, in general, the critical solution would also change. 

In the following section we discuss results from detailed studies of black hole threshold solutions 
generated from several distinct initial boson star states. Table |4~4~T1 summarizes the values of 4>o(0) 
that were used, the corresponding values of A3 required to generate a critical solution, and the 
figures that display results associated with the respective calculations. Since we will not dwell on 
this point below, we note that all of our calculations confirm the basic picture previously reported 
that the black holes that form just above threshold in this type of collapse generically have finite 
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mass (i.e. that the critical transition is Type I). 
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Table 4.1: Summary of parameters used to generate the results displayed in Figs. 14.414. 131 Listed 
are the figure number, central amplitude of the complex field, <^o(0), and the overall 
massless scalar amplitude factor, A\ fsee !4.10~|l . that generates a marginally-critical so- 
lution in each case. Other parameters defining the massless scalar initial profile (|4.109|) 
are held fixed at tq = 40, a = 5 for all simulations. 



4.4.2 Critical Phenomena: Results 

We start by examining results from a critically perturbed boson star having an unperturbed central 
field value </>o(0) = 0.05. As just described, the critical massless amplitude factor, A3 ~ 0.0032 was 
determined by performing a bisection search on ^3, to roughly machine precision. (Recall that 
each iteration in this search involves the solution of the time-dependent PDEs for the model for a 
specific value of A3, with all other parameters held fixed, and the criterion by which we adjust the 
bisection bracket is whether or not the simulation results in black hole formation. ) 

A series of snapshots of dM(t, r) jdr (where M(t, r) is the mass aspect function) for a marginally 
subcritical evolution is shown in Fig. 14.41 Full analysis of the results of this simulation indicate 
that the boson star enters what we identify as the critical state at t as 130, and remains in that 
state until t as 510. It is worth noting that the boson star actually completes its collapse into a 
more compact configuration well after the real scalar field has dispersed from the boson star region. 
We also note that the amount of time, r, spent in the critical state — r as 380 in this case — is a 
function of how closely the control parameter has been tuned to criticality. Specifically, we expect 
t to be linear in ln|Aa — A|| (see (|4.108t V and we will display evidence for this type of scaling 
below. 

Fig. l4.5l shows the time evolution of the central modulus of the complex scalar field for marginally 
subcritical evolutions generated from boson star initial states with (f>o(0) = 0.035,0.04 and 0.05. 
From the figure we can see that in all three cases the perturbed stars enter an excited, critical 
state at t as 100 and remain in that state for a finite time which is a function of </>o(0) (i.e. of the 
initial state). Additionally, at least for the cases <fo(0) = 0.035, cfo(0) = 0.04, the figure provides 
evidence that following the critical evolution phase, the excited stars relax to states characterized 
by large amplitude oscillations of the complex field. This behaviour will be examined in more detail 
below. Finally, also apparent in the plot are the smaller-amplitude oscillations during the periods 
of critical evolution. Previous work ~3 |~] indicated that these oscillations can be interpreted as 
excitations of the (stable) first harmonic mode of the unstable boson star that is acting as the 
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Figure 4.4: Critical evolution of a perturbed boson star with 4>o(0) = 0.05 and mass Mc = 
0.62Mp ; /m. This figure shows the time development of contributions to dM jdr from 
the complex (solid line) and real (dashed line) scalar fields. Note that the temporal 
spacing between successive snapshots is not constant — the time instants displayed 
have been chosen to illustrate the key features of the near-critical evolution. Also 
note that we have multiplied the value of dM/ dr for the real scalar field by a factor 
of 8 to aid in the visualization of that field's dynamics. The evolution begins with 
a stable boson star centered at the origin, and an in-going gaussian pulse (shell) of 
massless, real scalar field that is used to perturb the star. The overall amplitude 
factor, A3, of the initial real scalar field profile (see (|4.109|l . is the control parameter 
for generating the one-parameter family of solutions that interpolates through the 
black hole threshold. For the calculation shown here, A3 has been tuned to a critical 
value A3 w 0.0032 via a bisection search (and with a fractional precision of w 10~ 15 ). 
The other parameters defining the gaussian initial profile of the massless fields are 
ro = 40 and a = 5. The snapshots show that the real scalar field enters the region 
containing the bulk of boson star at t w 22, implodes through the origin at t w 45, 
leaves the boson star region at t « 70, and, finally, completely disperses from the 
computational domain at t m 100. The boson star enters the critical state at roughly 
the same time that the real field leaves the domain, and remains in that state for a 
period of time which is long compared to the crossing time of the massless field. At 
t w 510, the boson star departs from the critical state. 
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Figure 4.5: Time evolution of central value of the modulus of scalar field for subcritical evo- 
lution of perturbed boson stars. The figure shows the time evolution of \<f>(t, 0)| 
for marginally subcritical evolutions generated from boson star initial states with 
4>o(0) — 0.035,0.04 and 0.05. See the text for a description of key features of this 
plot. 



critical solution — the unstable fundamental mode is the one that determines whether or not the 
configuration will evolve to a black hole. Although we have not studied this matter in any detail, 
we assume that the same picture holds for our current calculations. 

The results from our simulations of critically perturbed boson stars are thus in agreement with 
previous studies |2j that identified the critical states as excited (primarily in the first harmonic 
mode), unstable boson stars. Following that work we can display an approximate correspondence 
between the initial boson stars and the critical solutions as show in Fig. 14. 61 The solid line shows the 
one-parameter family of static boson stars (parameterized as usual by 4>q(0)), where we have defined 
the radius, R, of a boson star so that M(R) = 0.99Madm = 0.99M(oc). The triangles indicate the 
initial stable boson star configurations, the squares indicate our best estimate of the corresponding 
unstable critical boson star states, and each arrow schematically depicts the transition between the 
two states that is induced by the perturbing scalar field. We note that to identify which unstable 



CHAPTER 4. BOSON STARS IN SPHERICAL SYMMETRY 



57 




5 10 15 20 

Radius R containing 0.99 M 

Figure 4.6: Transition of perturbed boson stars in critical evolutions. The solid curve shows the 
parametric mass vs radius plot of static boson stars (curve parameter, 0o(O)), where 
we have defined the stellar radius, R, so that M(R) = 0.99 M(oo) = 0.99 Madm- 
Triangles label the initial configurations, squares show the corresponding critical 
solutions (identified as one-mode-unstable boson stars with oscillations — largely in 
fundamental mode), and the dashed arrows represent schematically the transition 
between the initial and critical states. See the text for more details. 
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boson star is acting as the critical solution — which is equivalent to identifying an effective value 
of </>o(0) — we time average the central modulus of the complex field, \<f>(t, 0)| during the period of 
critical evolution. In addition, in accord with previous results, we observe that in all cases, the 
mass of the unstable critical state is larger than that of the progenitor boson star, indicating that 
a significant amount of mass-energy is extracted from the massless scalar field through its purely 
gravitational interaction with the complex field. 

As discussed previously, for both subcritical and supercritical simulations, the closer one tunes 
A3 to the critical value A3, the longer the perturbed star will persist in the critical state. Specifically, 
we observe scaling of the lifetime, r, of the critical evolution of the form 

r(A 3 )~- 7 hi|A3-A*|, (4.111) 

where we define the lifetime to be the lapse of coordinate time from the start of the evolution, 
t = 0, to the time of first detection of an apparent horizon, and where 7 is a scaling exponent that 
depends on which of the infinitely many one-mode unstable boson stars acts as the critical solution 
in the particular scenario being simulated. We note that the details of the definition of t are not 
important to the determination of 7 in (|4. 111(1 since 7 actually measures the differential in lifetime 
with respect to changes in A3 — A3, and this differential is insensitive to precisely how we define 
r, at least as A3 — > A3. In addition, we note that in using coordinate time in our definition of 
the scaling relationship (|4.111|l . we are defining the scaling with respect to proper time at spatial 
infinity. Another choice — arguably more natural — would be to define t in terms of the proper time 
measured by an observer at rest at r = (central proper time). Since the critical solutions are 
nearly static, the relation between these two different definitions of time would be a specific factor 
for each distinct value of ^o(O), and would thus lead to a 0o(O)-dependent "renormalization" of the 
scaling exponents, 7. 

Fig. 14.71 shows measured scaling laws from supercritical evolutions of perturbed boson stars 
defined by </>o(0) = 0.02,0.035,0.04 and 0.05. It is clear from these plots that, at least as A 3 — > 
A3, we have lifetime scaling of the form 1(4.111(1 . Estimated values of 7 — computed from linear 
least-squares fits to the plotted data— are 7 = 8.1,11,14,17 for O (O) = 0.02,0.035,0.04,0.05, 
respectively. We note that according to the now standard picture of critical collapse (see for 
example 29 ), each value of 7 can be identified with the reciprocal Lyapunov exponent (i.e. growth 
factors) of the single unstable mode associated with the corresponding critical solution. Again, 
the reason that we observe different values of 7 for different choices of initial boson star (different 
values of </>o(0)) is that distinct critical solutions are being generated in the various cases. That is, 
we cannot expect universality (with respect to initial data) in this case because the model admits 
an entire family of one-mode unstable solutions that sit at the threshold of black hole formation. 

4.4.3 Final Fate of Subcritical Evolutions 

In previous work on the problem of critically perturbed spherically symmetric boson stars ^ |2] , 
it was conjectured that the end state of subcritical evolution was characterized by dispersal of the 
boson star to large distances (relative to the size of the initial, stable star). This conjecture was 
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Figure 4.7: Measured lifetime scaling laws for critically perturbed boson stars. This figure shows 
the measured lifetimes of various near-critical evolutions of perturbed boson stars as 
a function of ln|A 3 - A*|, for cases with </>(0) = 0.02,0.035,0.04 and 0.05. Quoted 
scaling exponents, 7 (see (|4.111ll are computed from linear least-squares fits to the 
data. The apparent convergence of the data for different <^o(0) as In | A3 — A3 1 — > 
is not significant, as it reflects calculations far from criticality i.e. far from the 
In I A3 — A3 1 — > —00 limit. See the text for additional details. 



CHAPTER 4. BOSON STARS IN SPHERICAL SYMMETRY 



60 



at least partially influenced by the behaviour observed, for example, in the collapse of a massless 
scalar field [7], where subcritical evolutions do involve complete dispersal of the field. However, 
another key reason for what we claim is a misidentification of the true subcritical end-state was 
that the simulations described in 012] simply were not carried out for sufficient coordinate time to 
see the long-time behaviour. Our current simulations strongly suggest that subcritical evolutions 
lead to a "relaxation" of the critically perturbed state to something that approximates a boson 
star (not necessarily the original star) undergoing large amplitude oscillations. As argued in the 
next sub-section, these oscillations can largely be identified with the fundamental perturbative 
mode associated with the final boson star state. The numerical evidence also suggests that these 
oscillating configurations eventually re- collapse in general; a "prompt" re-collapse can be seen in 
the 0o (0) = 0.05 data in Fig.IPl 

Fig. 14.81 shows the long-time behaviour of max r (2M(t, r)/rs), \<p(t,0)\ and i/j(t,0) for a near- 
critically perturbed boson star (<j>o(0) = 0.04, ^ as 0.00342), where r$ = ip 2 r is the areal radius, 
and where, loosely speaking, 2M(t,r)/rs = 1 signals the surface of a black hole. Note that this 
is a subcritical evolution, so that a black hole does not form. As shown in more detail in previous 
figures, the boson star enters a critical state (well approximated by an unstable boson star) shortly 
after the real scalar field leaves the computational domain (t w 100), and while in the critical state, 
it oscillates with the frequency of the fundamental mode as computed from perturbation theory 
using the unstable boson star state as the background (see E]). At t w 370 the star leaves 
the more compact critical configuration, decreases in central density, expands in size, and starts to 
pulsate with a different frequency. Although at late time the oscillation amplitudes are much larger 
than those seen in the critical phase of evolution, we will show in the following section that the 
oscillations can nonetheless be largely attributed to excitations of the fundamental perturbative 
mode associated with the final boson star state. 

4.4.4 Perturbation Analysis of Subcritical Oscillations 

We now proceed to an application of perturbation theory to the oscillations seen in long-time 
evolutions of marginally subcritical configurations, such as those shown in Fig. 14.81 Here we 
follow \Tj\ an d d< and refer the interested readers to those sources for details of the approach 
that are not included here. In particular, we emphasize that we have not carried out the complete 
perturbation analysis ourselves, but are simply using a computer code provided by Hawley [Mj 
to analyze our current simulations. Nonetheless, to make contact between the perturbative and 
simulation results, it is useful to briefly review the setup of the perturbative problem. 

To formulate the equations for the perturbation analysis, we first rewrite the complex scalar 
field as: 

0(t,r) = (Vi(t,»0 +iMt,r))e~ luJt , (4.112) 

(Note that this representation is distinct from <f> = <f>± + i</> 2 , and the reader should be careful not 
to confuse the ^>'s used here with the conformal metric variable, ip.) Additionally, the spacctime 
metric is written in Schwarzschild-like (polar-areal) coordinates: 



ds 2 = -e^dt 2 + e x ^dr 2 + r 2 (d8 2 + sin 2 0dtp 2 ) 



(4.113) 
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Figure 4.8: Long time behaviour of subcritical evolution for <fi(0, 0) = 0.04. This figure shows the 
long-time behaviour of max r (2M(t,r)/rs), \4>{t, 0)| and ip(t,0) for a near-critically 
perturbed boson star (<fio(0) — 0.04, A3 « 0.00342), where r$ = ip 2 r is the areal 
radius, and where, loosely speaking, 2M(t,r)/rs = 1 signals the surface of a black 
hole. Note that this is a subcritical evolution, so that a black hole does not form. The 
figure provides evidence that the final state of subcritical evolution is characterized by 
large amplitude oscillations about something approximating a boson star on the stable 
branch, rather than dispersal of the complex field as suggested in Detailed 
calculation (see Sec. 14.4. 411 shows that the pulsation frequency is approximately the 
fundamental mode frequency computed from perturbation theory about a background 
stable boson star solution with <^>o(0) = 0.051. Also note the overall lower-frequency 
modulation of the post-critical-phase oscillations. This effect is not yet understood, 
although one possible explanation — namely that the envelope modulation represents 
"beating" of the fundamental and first harmonic modes — appears to be ruled out. 
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We further introduce four perturbation fields 6X(t, r),5v(t, r), 6tp\ (t, r) and 8ip2(t, r) which represent 
the perturbations about the equilibrium values Ao(r), v (r) , O (r) : 



X(t,r) = X (r)+5X(t,r), 

v{t,r) = v (r) + 5v{t,r) , 

^i(t,r) = 0o(r)(l + <tyi(t,r)) , 

V>2(£,r) = <l>o(r)6il>2(t,r) . 



(4.114) 
(4.115) 
(4.116) 
(4.117) 



With the above definitions we can write the coupled Einstein-Klein-Gordon field equations as a set 
of PDEs for the functions 5X, 5v, Sipi and 5ip2- With some manipulation we can then eliminate 5v 
and 5ip2 to produce a system of two coupled second-order PDEs for Sipi and SX: 
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(4.118) 
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(4.119) 



Note that these equations involve only second time derivatives (i.e. there are no terms involving 
Sipi or SX), and that they are linear in the second time derivatives. If we thus assume a harmonic 
time-dependence for the perturbed fields: 



<tyi(t,r) = 5Mr)e lat 
*Ai(t,r) = 5A!(r)e iCT * : 



(4.120) 
(4.121) 



then the equations for the perturbations contain a only in the form a 2 , and the sign of a 2 , as 
computed by solving a particular mode equation, determines the stability of that mode. (Note 
that the system can be shown to be sclf-adjoint so that the values of a 2 must be real.) If any of the 
values of a 2 are found to be negative, then the associated perturbations will grow and the boson 
star will be unstable. Moreover, as the eigenvalues form an infinite discrete ordered sequence, 
examining the fundamental radial mode tr 2 determines the overall stability of any particular star 
with respect to radial perturbations. 

In order to compare the simulation results with those given by perturbation theory, we first 
observe that there is a difference in the choice of the time coordinates used in the two calculations. 
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Specifically, in the perturbative analysis ^| |2] , the lapse was chosen to be unity at the origin, so 
we have 



T 2 

2 



pcrturbativc (X 



simulation 



We also note that there is a factor of 2 difference in the definitions of T^ v used in the two calcula- 
tions, and that the definition of the complex field, (f>(t,r), in the perturbative calculation includes 
a factor of y/8n. We thus have 



— > V47T(/> 

perturbative 



simulation 



The numerical technique for obtaining the fundamental mode and first harmonic mode frequencies 
of boson stars has already been described in [2] and will not be repeated here; again we will 
simply quote and use results from that study. From Fig. 14. 81 we note that there are 20 oscillations 
between t = 620 and t = 5200, giving a period T w 229. Hence we have an oscillation frequency 
fj = 2ir/T w 0.0274. The time average of the lapse function (a(t, 0)) in the interval is 0.78, and 
so a 1 jo? s» 0.0013. We also compute the time average of 4>(t,0) in the interval, and use the 
resulting value to identify the stable boson star solution about which we perform the perturbation 
analysis. We find (<fo(t,0)) ~ °- 051 x V4vr = 0.18. For a boson star with <j) (0) = 0.18, the 
perturbative calculations (see Fig. 7 of predict a 2 = 0.0014, which is in reasonable agreement 
with the simulation results. Hence the oscillations that occur in the post-critical regime appear to 
be largely fundamental mode oscillations of a final-state, stable, boson star. We also remark that 
since the oscillations are of such large amplitude, it does not appear possible to precisely identify 
an effective background state (i.e. an effective value of </>o(0)), so the level of agreement in the 
oscillation frequencies is possibly as good as one could expect. 

We have also attempted to understand the nature of the slower frequency oscillation which mod- 
ulates the envelope of the fundamental-mode oscillations, as is visible in Fig. 14.81 One hypothesis 
is that this modulation represents a "beating" effect of two oscillatory modes with frequencies u\ 
and o 2 such that \o\ — 02 1 <C cri,2- We observe that period of the envelope of oscillation is roughly 
T « 2000, corresponding to a squared-frequency a — 2tt/T w 0.00314. The average of the lapse 
function during one cycle is roughly 0.78, so a 2 /a 2 w 1.62 x 10~ 5 . However, the first harmonic 
frequency from perturbation analysis for (4>o{t, 0)) » 0.18 is a\ = 0.0105, much different from the 
fundamental mode a 2 = 0.0014. Therefore the beating cannot be produced from the superposition 
of the fundamental and first harmonic modes. 



4.5 Black Hole Excision for Boson Stars 

In the last section of this chapter we describe our application of black hole excision techniques in 
the context of gravitational collapse of boson stars. As described in Chap. |3| black hole excision 
methods, whereby the interiors of black holes are excluded from the computational domain, have 
been the subject of intense study in numerical relativity over the past decade. We view our current 
implementation as a further proof-of-principle for the strategy, as well as providing some indications 
of the subtle interactions between inner (excision) and outer boundaries that can arise when using 



CHAPTER 4. BOSON STARS IN SPHERICAL SYMMETRY 



64 



the excision method. We also emphasize that we have not yet used calculations with excision in 
order to generate new results of physical significance. 

To study excision techniques, we naturally need to evolve spacetimes which eventually develop 
black holes. Conveniently, we can consider the same model and initial data configurations used in 
our previously-described study of critical phenomena, provided that a sufficiently large amplitude 
factor, A3, for the scalar field is used. In particular, we need to choose A 3 large enough so that 
the black hole that forms is large compared to the spatial discretization scale, Ar, used in the 
simulation. For the purposes of illustrating our implementation of excision, we consider a boson 
star with </>o(0) = 0.01, and perturb it with an imploding massless scalar field having an initial 
profile of the form (|4.1U9|I with A3 — 0.041, ro = 40 and a — 5. In fact, A3 is large enough in this 
instance that almost all of the mass associated with both the real and complex fields ends up in 
the final black hole. 

Fig. l4.9l shows a series of snapshots of 2M(t, r)/rg during the resulting evolution, where M(t, r) 
is the mass aspect function and rg — tp 2 r is the areal coordinate (r is the isotropic radial coordi- 
nate). During the evolution we monitor (|B.22|) to detect the appearance of an apparent horizon — in 
this case an apparent horizon is first detected at t — 80.51 and at r = 9.67. Assuming cosmic cen- 
sorship, the region r < 9.67 lies within a black hole for all future times so can safely be excised 
from the simulation domain for t > 80.51. (In Fig. 14.91 the region being excised lies to the left of 
the dashed vertical line in each sub-plot.) Fig. 14.101 shows the relative time rate of change of the 
mass aspect function at outer boundary, c> t M(i, r max )/Af (i, r max ), for the same evolution. At late 
time (t > 100) the relative change in the mass aspect function at outer boundary is of the order of 
10~ 6 , indicating almost no mass is passing through the outer boundary. 

Once excision is enabled, the origin is no longer part of the computational domain and thus 
the inner boundary (regularity) conditions previously imposed at r = for the variables a, tp and 
K r r , must be replaced by conditions imposed at the excision surface. In the case of the dynamical 
geometrical variables, tp and K r r , we can use discrete versions of the evolution equations l|B.23(l 
and (|B.24|I to update values for tp and K r r on the excision boundary. For the case of the lapse 
function, a, we have found that a simple strategy, whereby we 'freeze" the inner value of the lapse 
to whatever its value is when excision starts, leads to stable, convergent results. (Also note that 
since our simulations do appear to converge using our simple strategy, we have not experimented 
with the condition on a(r e , t), where r e is the r-coordinate of the excision surface. So, for example, 
we have not investigated whether the condition on a is perturbatively stable.) Finally, although 
the shift component, (3, also satisfies a regularity condition, the boundary condition for l|4.49|l . 
which governs the shift, is naturally set at large r and thus we do not need to specify a condition 
on (3 at the excision surface. 

With regards to our simulations with excision, a key observation is that after the coupled boson 
star/real scalar field configuration promptly collapses to form a black hole, the region exterior to 
the black hole quickly settles down to an essentially Schwarzschild spacetime. In particular, the 
mass aspect function outside the black hole approaches an almost-constant function, and a stable 
simulation of the black hole spacetime can be carried out for very long integration times. 

However promising the excision technique appears, careful examination of the mass aspect 
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Figure 4.9: Evolution of a highly perturbed boson star (4>o(0) = 0.01) using black hole excision. 

This figure shows snapshots of 2M(t, r)/rs, where r$ = ijj 2 r is the areal radius, and 
where, loosely speaking, 2M(t,r)/rs — 1 signals the surface of a black hole. In 
snapshots with a vertical dashed line, the region from r = to the radial position of 
the dashed line is excised from the computational domain. The evolution starts with 
a stable boson star centered at the origin, and a high-amplitude, in-going pulse of real 
scalar field of the form l)4.10QJl (A3 = 0.041, a = 5 and r = 40). At t w 80 the real field 
reaches the origin, and an apparent horizon containing most of the combined mass 
of the real and complex scalar fields forms at t = 80.51. At late times, the exterior 
solution settles down to a configuration very close to a static Schwarzschild spacetime. 
Note that the temporal spacing between successive snapshots is not constant-the time 
instants displayed have been chosen to illustrate the key features of the simulation. 



CHAPTER 4. BOSON STARS IN SPHERICAL SYMMETRY 



6(5 



5 -0.0002 



-0.0004 



S -0.0006 



-0.0008 



-0.001 



100 



200 



300 



Figure 4.10: Relative time rate of change of the mass aspect function at outer boundary, 
dtM(t, r max )/M(t, r m ax), for the evolution of a highly perturbed boson star (</>o(0) = 
0.01) using black hole excision. The figure shows the result for the same evolution in 
Fig. 14.91 with r max = 150 and Ar w 0.146. At late time t > 100 the relative change 
in the mass aspect function at outer boundary is of the order of 10~ 6 , indicating 
almost no mass is passing through the outer boundary. 
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function over sufficiently long integration times does reveal certain numerical artifacts. For example, 
Fig. 14. 1 II plots the mass aspect function at the outer boundary of the computational domain, 
r = r max = 150, against time, for < t < 360. During any period where there is no flux of real or 
complex mass/energy crossing the outer boundary, we should have M(i, r max ) = const. From the 
figure we see that M(t, r max ) decreases initially, mainly due to the escape of real scalar field from 
the computational domain. At t w 80, following the first detection of an apparent horizon, and the 
enabling of excision, there is a relatively sharp, and transient decrease in M(t, r max ). However, as 
can also be seen in the figure, all indications are that this particular decrease "converges away" 
in the continuum limit Ar — > 0. Following the transient phase, although the mass generically 
remains roughly constant, there is always a small rate of increase in M(i,r max ) which is largely 
independent of the grid resolution, Ar. This is surely an artifact — since no mass is physically added 
to the system, we claim that the observed mass increase must be a purely numerical/computational 
effect, and we have thus investigated possible remedies. 

Because our simulations are performed on a finite radial domain, < r < r max , it is nat- 
ural to consider the effect of the placement of the outer boundary of the domain, r max , on the 
results. Fig. 14.121 plots M(t, r max ) vs t from simulations using the same initial data described 
in Fig. 14.111 but where r max has been doubled to r max = 300. From this figure it can be 
seen that the transient mass decrease immediately following the enabling of excision still con- 
verges away as Ar — * 0. In addition, it is apparent that the level of mass increase following 
the transient phase is substantially reduced from the r max = 150 calculations. To illustrate this 
point more clearly, Fig. 14.131 shows the deviation of the mass aspect function M(t, r max ) from 
Af(190,r max ) vs time t for r max = 150,300 and 600. The deviation is measured relative to the 
mass aspect function at t — 190, by which time the transients associated with the enabling of 
excision have settled down. The figure shows that the deviation is smaller when we increase 
the location of the outer boundary from ?* max = 150 to r max = 300 and r max = 600. The calcu- 
lated value for (AM (t, 150) - AM (t, 300)) / (AM(i, 300) - AM(i, 600)) « 5, where AM (t, r max ) = 
M(t, r max ) — M(190, r max ), which indicates the rate of convergence is roughly r~ ax . 

We thus conjecture that the most dominant contribution to the small amount of post-excision 
"mass inflation" we observe is due to our (approximate) asymptotic treatment of the boundary 
conditions. We have done some preliminary investigations of improvements to the outer boundary 
conditions but have no significant findings to report at this time. 
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Figure 4.11: The mass aspect function M(i,r max ) vs t for r max = 150 from simulations with 
excision, and computed using three different finite difference resolutions, Ar, in a 
4:2:1 ratio. For resolution Ar = 0.146, (the coarse calculation), the mass decreases 
quite suddenly (and by a few percent) at t ss 80, at which time a black hole is de- 
tected and a region of spacetime is excised. Following this early transient associated 
with the enabling of excision, the mass remains roughly constant (relative to the 
transient fluctuation), although a slight constant rate of mass increase is observed. 
As resolution increases (Ar = 0.073, Ar = 0.037), the early transient mass decrease 
apparently converges away (possibly only as Ar). However, the late-time mass in- 
crease, does not seem to vanish as Ar — > as can be seen from the fact that the 
slope of the plots for t > 100 remains roughly constant as the resolution is varied. 
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Figure 4.12: The mass aspect function M(t,r max ) vs t for r max = 300 from simulations with 
excision, and computed using three different finite difference resolutions, Ar, in a 
4:2:1 ratio. As in Fig. 14. Ill with increasing resolution, the magnitude of the early 
(t ~ 80) transient mass drop apparently converges away. Moreover, the mass is 
more nearly conserved at late times in the calculation, relative to the simulation 
shown in Fig. 14.111 In other words, the spurious mass increase seems to be due to 
our treatment of the outer boundary conditions, and we can apparently suppress 
the resulting numerical artifacts by increasing the size of the computational domain. 
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Figure 4.13: The deviation of the mass aspect function M(t, r max ) from M(190, r max ) vs t for 
'"max = 150, 300 and 600 from simulations with excision, and computed using three 
different finite difference resolutions, Ar, in a 4:2:1 ratio. The deviation is measured 
relative to the mass aspect function at t = 190, by which time the transients asso- 
ciated with the enabling of excision have settled down (see Fig. 14. Ill and Fig. l4.12|) . 
The deviation is smaller when we increase the location of the outer boundary from 
''max = 150 to r max = 300 and 7- max = 600. The rate of convergence is roughly r"^. 
See the text for further details. 



CHAPTER 5. BOSON STARS IN AXISYMMETRY 



71 



CHAPTER 5 

BOSON STARS IN AXISYMMETRY 

We now consider the evolution of the coupled Einstein-Klein-Gordon field equations in axisymme- 
try. As time-dependent, spherically symmetric spacetimc is the most natural 1 space + 1 time- 
dimensional arena in which to approximate star-like objects, so is time-dependent axisymmctric 
spacetime the most natural 2+1 context. Crucially, the generalization to axisymmetry allows us 
to consider time-independent configurations with net angular momentum, which, for the case of 
boson stars turn out to be surprisingly non-trivial to construct. In addition, with an axisymmetric 
code we can simulate scenarios such as the head-on collision of two or more boson stars, as well as 
axisymmetric (i.e. non-spherically symmetric) critical collapse. This chapter concerns the construc- 
tion of stationary (time-independent 1 ), rotating boson star solutions (which requires the solution 
of an elliptic quadratic eigenvalue problem in two spatial dimensions), as well as the construction 
of time-dependent axisymmetric solutions representing the evolution of generic (in principle) ax- 
isymmetric initial data for a complex scalar field. Unfortunately, the stationary solutions that we 
construct in the first case are not time-independent in the coordinate system used in the latter 2 
and we are thus unable to directly study the evolution of rotating boson stars by combining the 
two calculations. We are, however, able to reproduce and extend previous calculations of rotating 
boson star configurations, as well as to study critical collapse of driven boson stars in axisymmetry 
through a completely original set of calculations. Some of the work presented in this chapter was 
done in collaboration with Choi |75| , who performed numerous tests on the solutions for stationary 
rotating boson stars in the Newtonian limit, as well as tests of the modified graxi code. Results 
of these tests are available on-line [JSj and will not be further described here. 

Before proceeding to a summary of the remainder of the chapter, we present an overview of some 
key features of general relativistic, rotating boson stars (and previous work on the subject), and 
then an overview of related work on axisymmetric critical collapse as well as our general approach 
to that problem. 

5.0.1 Rotating Boson Stars 

The question of existence of rotating fluid stars (solutions of the stationary, axisymmetric Ein- 
stein/hydrodynamics equations) has a clear answer on physical grounds. The observation of ro- 

1 We emphasize that when we speak of time-independent solutions in this chapter, we mean that the geometry is 
time- independent. In general, the complex scalar field, rp(t, x), will always be time-dependent, even if that dependence 
is "trivially" harmonic. Similarly, when we consider rotating boson stars, and speak of axially symmetric solutions, 
we will also mean that the geometry is axially-symmetric; again, the particular ansatz that we will adopt for a 
rotating boson star will result in explicit ip dependence in cf>(t, r, 6, ip). 

2 Specifically, the coordinate system (i,p~,z,(p) currently used in graxi for evolving complex scalar fields with 
rotation 1521 is not adapted to the timelike Killing vector of our stationary solutions, and so a non-trivial coordinate 
transformation — itself involving solution of a complex set of non-linear PDEs — must be performed to determine the 
initial data [<f>(t, p, z, ip), il(t, p, z, p)\ required by graxi. 



CHAPTER 5. BOSON STARS IN AXISYMMETRY 



72 



tating neutron stars — for example, in the form of pulsars — tells us that relativistic, rotating fluid 
stars naturally occur in our universe, and the theoretical study of such objects has been an active 
area of research for decades (for example, see [221 E3 El HH1 EH EDI) ■ Interestingly, it was not 
until fairly recently that the theoretical existence of relativistic, rotating boson stars was demon- 
strated [221 EH This is perhaps not too surprising when one considers the lack of astrophysical 
importance (to date at least!) of boson stars. In addition, at a more technical level, the regularity 
condition for a rotating bosonic star is different from that for a rotating fcrmionic star, a fact which 
leads to completely different topological shapes for the rotating stationary states in the two cases. 
For example, as has already been mentioned, for boson stars with non-zero angular momentum, the 
stationary solutions have mass-energy concentrated in a roughly toroidal (rather than spheroidal) 
configuration. 

Schunck & Mielke , as well as Yoshida & Eriguchi [2Hj have previously constructed general 
relativistic, stationary, rotating boson star solutions using numerical means. One of the more 
interesting features of such solutions is that angular momentum is quantized (or perhaps more 
properly, discretized). Specifically, to generate time-independent, rotating configurations, one is 
naturally lead to an ansatz (in spherical-polar coordinates) 3 

W,r,0,<p) — Mr,0)e iut+ik,p • (5.1) 

Now, u> can in principle take on arbitrary real values; k, however, must be integral for (f> to be single 
valued. Since the specific choice of k affects boundary conditions at r = 0, the ansatz (|5.1|) itself 
naturally leads us to consider a series of (elliptic eigenvalue) problems for k = 1,2, • •• (the case 
k = describes the spherical boson stars of Chap.0J). 

As in the spherical case, k = 0, for any value of k = 1, 2, • • • , u) will be the eigenvalue w(</>/ fc )), 
where is a generalization of the central modulus value, <po, conventionally used to parametrize 
the family of stars in spherical symmetry. In other words, as was the case in spherical symmetry, 
for each value of the angular momentum parameter (or azimuthal quantum number), k, we expect 
to find an entire family of rotating stars, which will be parametrized by 

Given the clear analogy between l|5.1|l and axisymmetric bound states in quantum mechanics, 
it is perhaps not surprising that rotating boson stars have distinctive topological characteristics. 
In fact, given the ansatz l|5.1|l one can show that regularity as r — ► implies 

lim Mr, 0; k) = r k f k (9) + 0(r k+2 ) . (5.2) 

Thus, for k > (i.e. the non-spherically symmetric cases) we have that <j> and its first k — 1 radial 
derivatives vanish on axis, and are relatively small close to the axis. This means that for k > 1, 
the boson stars might be more accurately termed boson tori. 

Solutions for rotating boson stars using the full general relativistic equations were first obtained 
by Schunck & Mielke [2(11 18 1). Specifically, configurations for 2 < k < 10 as well as for k — 500 were 
found. However, as pointed out by Yoshida & Eriguchi [2Hj the solutions obtained were all close to 

3 As noted in the introduction, this ansatz was apparently first written down, for the case of Newtonian boson 
stars, by Silveira & de Sousa 1241 . and is not the most general ansatz that could be considered. 
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Newtonian; for instance, in [SI] a k = 1 solution is presented which has M = 0.07331Mpj/m, a mass 
that is only about 6% of the maximum possible mass for the k = 1 family. Overall we estimate that 
none of the solutions constructed by Schunck & Mielke (for any of the values of k considered) have 
masses in excess of 10% of the maximum mass of the corresponding family. As another example 
of the weak-field nature of the Schunck & Mielke solutions, the quasi-Newtonian potentials that 
are displayed have almost no angular dependence, whereas the potentials corresponding to highly 
relativistic stars do. Moreover, Schunck & Mielke determined only a few configurations for any 
value of k, and thus did not study any specific family of solutions in sufficient detail to determine, 
for example, whether a "Chandrasekhar limit" existed, as is the case for fluid stars, as well as for 
spherically symmetric boson stars. Finally, some of the solution plots in (such as Figs. 6.4 and 
6.6) are not very smooth, indicating that the calculations may have been underresolved. 

The work by Yoshida & Eriguchi [2H| represented a significant advance over Schunck & Mielke's 
in several respects. The solutions they constructed were smooth, and they were able to compute 
the entire family of solutions for the case k = 1. Further, they showed that the plot of ^adm(^Cc)) 
(where is again the family parameter), exhibited a maximum, and concluded that the maximum 
mass for a k = 1 rotating boson star is M max = l.SlM^/m. 

Despite the considerable success of the calculations described in , we have identified several 
unsatisfactory points about that work as well, and have at least partially ameliorated the identified 
weaknesses in our current research. First, the code used in |28| breaks down before the maximum- 
mass star for k = 2 can be computed, whereas the method we describe below can be used to compute 
the complete family for the k — 2 case. Second, Yoshida & Eriguchi's numerical approach uses a 
so-called "self-consistent field" method in which the solutions are represented as two-dimensional 
surface integrals at every mesh point (the mesh itself is two-dimensional); this means the solution 
algorithm is quite expensive. In contrast, we adopt a finite-difference approach that uses homotopic 
and multigrid techniques to solve the resulting non-linear, elliptic eigenvalue problem, and although 
we have no way of performing direct timing estimates, we feel that our approach is likely to be 
significantly more efficient computationally. 

Finally, although our calculations are currently limited by regularity problems (so that, for 
example, we cannot compute solutions for k = 3, 4, • • • ), we have reformulated our basic equations 
in a manner that we feel will largely solve those problems (see Sec. 15.2.61 . Although work on the 
modification of our algorithm to solve this new system is not complete, we feel that with it we 
will be able to compute families of solutions for k > 2, and, importantly, will be able to do proper 
convergence tests to assess the accuracy of the results. The lack of convergence tests or other means 
of estimating the accuracy of solutions is another shortcoming of both the Schunck & Mielke and 
Yoshida & Eriguchi work, and, indeed, our regularity problems inhibit the complete convergence 
testing of our current results as well. We note that regularity problems near coordinate singularities 
in curvilinear coordinates have often been encountered in calculation in numerical relativity. In this 
case, we have no rigorous understanding of the origin of these problems; however, we have strong 
evidence that the solutions of the discrete equations are not regular as r — > (i.e. our solution of 
the non- linear eigenvalue problem is currently unstable). 



CHAPTER 5. BOSON STARS IN AXISYMMETRY 



74 



5.0.2 Axisymmetric Critical Phenomena 

An investigation of axisymmetric critical phenomena very similar to the one described below was 
carried out by Rousseau |31|. who studied critical phenomena associated with driven collapse of 
spherically symmetric boson stars in the context of the so-called conformally flat (or Isenberg- 
Mathews- Wilson) approximation to general relativity. Using a massless scalar field as the per- 
turbing agent, Rousseau was able to verify the existence of families of solutions that interpolate 
across the black hole threshold (the conformally flat approximation admits solutions with apparent 
horizons), and found some evidence of Type I behaviour — with the critical solutions being unstable 
boson stars — and indications of lifetime scaling of the form l]5.45[l . However, Rousseau's calcula- 
tions were not fully relativistic, and were severely resolution-limited, so the results are neither very 
complete nor accurate. 

All of the dynamical studies described below have been performed using a code, graxi, pre- 
viously developed by Choptuik, Hirschmann Liebling & Pretorius [SJ |§] , that solves the coupled 
Einstein-Klein-Gordon equations in axisymmetry. I was responsible for the implementation of the 
initial data interface for boson stars, as well as for the execution and analysis of the numerical 
experiments reported in this Chapter. We note that the dynamical calculations that we describe in 
Sec. 15.31 are quite well resolved and accurate (estimated local accuracy of a few percent) by virtue 
of the fact that graxi incorporates Berger-and-Oliger-style adaptive mesh refinement. 

We have used graxi to investigate Type I critical behaviour of axisymmetric boson stars using 
two novel strategies. The first involves the head-on collision of two identical spherically symmetric 
boson stars, where the family parameter is the initial momentum, p z , imparted to each star. The 
second involves perturbation of a single spherically symmetric boson star via the implosion of a 
non-spherically-symmctric pulse of massless scalar radiation. In this case the family parameter is 
the overall amplitude of the initial massless scalar pulse. In both instances we find strong evidence 
for Type I behaviour, as well as strong evidence for lifetime scaling of the form 15.45f) . In addition, 
in some of the collision simulations we observe interesting "solitonic" behaviour that has previously 
been seen in the collision of Newtonian boson stars [12] . 

5.0.3 Outline 

The remainder of this chapter is organized as follows. We start in Sec. 15.11 with a review of the 
most general geometrical notion of "time-independence" for a generic spacetime (stationarity) and 
a discussion of a suitably general, symmetry-adapted coordinate system in which to study axisym- 
metric, stationary solutions of the coupled Einstein-Klein-Gordon equations. We then proceed in 
Sec. 15.21 to a discussion of an approach to the construction of solutions representing general rel- 
ativistic stationary boson stars with angular momentum. As with the determination of spherical 
boson star data per se, this involves the solution of the initial value (constraint) equations for the 
gravitational field, but is significantly complicated by the fact that we must simultaneously solve 
a non-linear elliptic eigenvalue problem (an elliptic PDE) for the matter field. In contrast to the 
previous work of Yoshida & Eriguchi , which employed a so-called "self-consistent field" tech- 
nique, we proceed via finite-difference solution of the stationary PDEs for rotating boson stars and 
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outline some of the computational challenges that result in Sec. 15.2.21 In Sec. l5.2.3l we then discuss 
details of the multigrid-based solver that we have designed and implemented to solve the eigenvalue 
problem. Results for rotating boson stars generated using the solver are presented in Sec. 15.2.41 
and we follow with remarks concerning the performance of our method and possible enhancements 
in Sec. l5~2~5l 

As already mentioned, our solver for rotating boson stars is currently limited by regularity 
problems, and we are, in fact, unable to compute solutions for k = 3,4, One promising 
approach to resolving this difficulty is to redefine variables (both scalar field and geometric) so that 
regularity in the discrete domain is more automatically enforced. To that end we present a candidate 
reformulation of the equations in Sec. 15.2.61 Although we have not yet implemented a solver for 
a discretization of this new system, we are confident that this task will be straightforward, given 
the computational techniques that we have already developed. We thus hope that the incomplete 
work described in Sec. 15. 2. 61 will soon give rise to a solver that will work for k > 2. 

In Sec. l5.3l we turn attention to the dynamics of boson stars. We discuss the results from a study 
of head-on collisions of two boson stars in Sec. 15. 3. 21 and finally, in Sec. 15.3.51 we summarize results 
from axisymmetric generalizations of the spherical critical collapse studies previously described in 
Sec. El 

5.1 Stationary, Axisymmetric Spacetime 

As mentioned above, the concept of stationarity in general relativity is meant to capture the 
most general geometric notion of "time-independence" . A spacetime is stationary if it admits a 
timelike Killing vector field, and operationally, we will equivalently define a stationary spacetime 
as one for which, at least locally, there exists a coordinate system in which all metric components 
are independent of the time coordinate, t. 4 In addition, we now wish to restrict attention to 
axisymmetric spacetimes, so that the spacetime also has a spacelike, azimuthal Killing vector (with 
closed orbits). Again, we will assume that this means that we can choose a coordinate system 
with one of the coordinates, tp, adapted to the symmetry, so that none of the metric coefficients in 
that coordinate system depend on (p. In other words, for stationary, axisymmetric spacetimes, if 
we adopt "natural" (i.e. symmetry-adapted) spherical-polar coordinates, (t,r,6,<p), the metric will 
depend only on r and 6: 

9»v(t,r,6,4>) — ► g^(r,0) . 

Since we are considering rotating configurations, our spacetimes will in general not exhibit time 
reversal symmetry as t — > — t will result in a configuration representing stars rotating in a sense 
opposite to the original. (In other words, stationary spacetimes are not necessarily static.) However, 
our spacetimes should be invariant if we simultaneously reverse both time and the azimuthal angle, 
i.e. the metric should have the discrete symmetry (t, ip) — ► (— t, —<p). In the general expression for 
the line-element in symmetry- adapted coordinates this fact immediately excludes the dtdr, dtd8, 

4 This will be true if the time coordinate is adapted to the timclikc Killing vector field \ a = (d/dt) a . In this case 
the Killing equation becomes C x g^ v = dg^u/dt = 0, i.e. the metric components are time-independent. 
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drdip and d6dip terms as they do not preserve this discrete symmetry. Finally, we can perform 
a coordinate transformation |S2] such that the metric of the 2-dimensional r-9 subspace has the 
diagonal form 

^ (dr 2 + r 2 d9 2 ) , 

as long as the subspace has a positive-definite (or negative-definite) signature. With these con- 
siderations we can now write a suitably general form for the metric of a stationary, axisymmetric 
spacetimc: 

ds 2 = -a 2 dt 2 + u (fidt + dip) 2 + ?/> 4 (dr 2 + r 2 d9 2 ) . (5.3) 

To summarize, by an appropriate choice of coordinates, we can reduce the number of non-trivial 
metric components that arc needed to describe stationary, axisymmetric spacetimes to four; here 
we identify those four quantities with the functions a,f3,ip and u. In what follows, we will also 
restrict attention to those stationary, axisymmetric configurations with equatorial symmetry (i.e. 
those that are symmetric under 9 — > tt — 9) . We thus can compute using the reduced angular range 
< 9 < it/2. 

5.2 The Initial Value Problem for Rotating Boson Stars 

5.2.1 The System of Equations in Quasi-Isotropic Coordinates (System 
A) 

Paralleling what was done in the case of spherical symmetry in the previous chapter, we have 
adopted what we will call the quasi-isotropic (spatial) coordinate condition in choosing the form (|5.3I) 
for the metric of a general axisymmetric, stationary spacetime. In fact, it is convenient to take u 
in H5.3fl to be of the form 

u = t/>V sin 2 9e 2a , (5.4) 

so that the deviation of a(r, 9) from zero reflects the deviation of the 3-metric from conformal 
flatness. Thus, the spacetime metric now becomes 

ds 2 = ~a 2 dt 2 + ^r 2 sin 2 9e 2a (J3dt + dip) 2 + t/> 4 (dr 2 + r 2 d9 2 ) , (5.5) 

or 

ds 2 = (-a 2 + ^r 2 sin 2 9e 2a /3 2 ) dt 2 + ^ (dr 2 + r 2 d9 2 + r 2 sin 2 9e 2 ° dp?) + 2t/>V sin 2 9e 2a f3dtd<p 

(5.6) 

where, again, a, f3, ip and a are functions of r, 9 only. To generate stationary solutions (having 
angular momentum in general) we adopt the following ansatz for (j)(t, r, 9, ip): 



^(t,r,9,ip) = Mr,0)e~ l(ut+kv) 



(5.7) 
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As discussed above, and is the case for a quantum mechanical wave function, k must take on integer 
values, since <p must be single valued with respect to (p: 



4(t, r, 9, if) = <f>(t, r,9,cp + 2tt) . 



(5.8) 



As we have also previously remarked, this "quantization" of angular momentum is one of the more 
novel features of rotating boson star solutions relative to rotating fluid configurations. 

We can now proceed to derive the PDEs for stationary boson stars in quasi-isotropic coordinates 
from the Hamiltonian and momentum constraints, the Klein-Gordon equation, the ansatz 1)5.70 . 
and the assumption of stationarity. Since the derivation is lengthy but straightforward, we will 
simply state the resulting equations. For numerical convenience we choose n — 1 instead of the 
usual convention k — 8ir (see (|G.l|l ) as in the previous chapter. 

First, from the Klein- Gordon equation we have the following 
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and where for simplicity of notation we have dropped the subscript "0" so that <f>o — > <f>. Note 
that from regularity considerations, the leading order behaviour of 4>Q(r,9) as r — > is cj>o(r,8) = 
r k (f>ok(9) + 0{r k+2 ). Therefore, in principle the equilibrium solutions can be parametrized by the 
central value of the kth radial derivative of the scalar field d k (j>o(Q,9)/dr k . In fact, numerically 
it may be more reasonable to define a new function <fr (r,9) such that (/>o(r,6) = r k <fr (r,9) (see 
Sec. l5.2.5)l . However, our current approach simply use 4>o(r, 9) as one of our fundamental variables. 
The Hamiltonian constraint gives 
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The r and 9 components of the momentum constraint vanish identically, while the ^-component 
of the momentum constraint gives 



P,rr-\ 5~+ - + 3(7. r h 6— P,r + 3cotfc> + 3(7,0 h 6— —5- = -167T0 ; „ 

r z \r a ^ J \ a ip J r z r 2 sin 9e 2a 

(5.11) 

From the assumption that the solutions are stationary, and because we work in a coordinate system 
adapted to the timelike Killing vector, the metric components must all be ^-independent. As a 
result, the left and right hand sides of the evolution equations for the metric l|2.17|l and extrinsic 
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curvature (|2.20|) must vanish. We can thus construct equations for the remaining variables, a 
and er, by taking linear combinations of the right hand sides of those equations. The combination 
K r r + K e e + K<p v = K = gives 
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Finally, the combination K r r + K 9 g — = gives 
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(5.13) 



Again, the leading order behaviour of cr(r, 0) as r — ► is cr(r, ff) — r 2 ai{&) + 0(r 4 ). (To properly 
impose the regularity condition it would be more appropriate to define a new function a(r, 9) such 
that a(r,9) = ra(r,6).) 

(|5.9|l - l|5.13() comprise our system of equations for stationary, axisymmetric boson stars. For 
any specific choice of k — 0, ±1, ±2, • • • these equations constitute a coupled, nonlinear eigenvalue 
problem for the five unknown functions a,j3,ip,a and <b (with eigenvalue w). For convenience we 
will refer to (|5.9|) - H5.13|) as System A. 



5.2.2 Challenges and Comments 

Before discussing our strategies for solving the above system, we would like to address some key 
challenges we face in solving System A. The discussion here motivates the particular computational 
strategies that we adopt in the following subsection. 

1. System A is nonlinear: numerical solution of nonlinear systems typically proceeds via iter- 
ation. Success of an iterative method is often contingent on the existence of a good initial 
guess with which to start the iteration. The definition of a "good initial guess" depends on 
both the method of solution being used, as well as the problem being solved. Our experience 
with the numerical solution of System A suggests that the basin of attraction for a given 
solution is often very small. Thus, we need to have the means of providing initial estimates 
which are already quite close to the solution. 

2. System A is a 2D eigenvalue problem: We could proceed to solve this problem using a direct 
method that generalizes the strategy outlined in Sec l4.3.3l for the case of spherically symmetric 
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boson stars. This would ultimately involve using linear algebra software to solve a general full 
matrix eigenvalue problem, which would require 0(N 3 ) operations, where N is the dimension 
of the matrix. For a 2D problem with modest resolution, say N r = Ng = 100, the dimension 
of the matrix is N — 10 4 , so the eigenvalue solution would require O(10 12 ) operations. This 
is obviously an expensive approach, even with current computers, not to mention that we 
would still have to iterate the eigenvalue solution. The direct technique also requires a large 
amount of memory to store the full matrix. We therefore need to focus on devising a more 
efficient algorithm for solving the eigenvalue problem. 

3. System A is a quadratic eigenvalue problem: the eigenvalue ui appears both linearly and 
quadratically in the equations of System A. Our solution algorithm must therefore be able 
to treat quadratic eigenvalue problems. 

4. System A has mixed (Robin) boundary conditions when the outer boundary of the compu- 
tational domain is located at a finite radius. Dirichlet boundary condition are the simplest 
conditions to implement, but we may need to impose conditions such as (|4.80|) . 

5.2.3 Numerical Strategies 

Having reviewed the difficulties we encountered in solving System A, we will now outline the 
solution strategies we have adopted. Since the specific form of the equations we eventually solve is 
lengthy, we focus here on a high-level description of the solution technique, and refer the interested 
reader to App.[0for more details. 

1. Compactification of coordinates: 

We introduce a compactified radial coordinate £, as well as a new angular coordinate, s: 

C = ~\~±L ' ( 5 - 14 ) 
1 + r 

s = cos9. (5.15) 

Thus we have £ <E [0, 1] and s 6 [0, 1]. This transformation simplifies some of the boundary 
conditions. For instance 14.80fl becomes 

VKC,s)lf=i = l- (5-16) 

Importantly, compactification not only allows us to work with (algebraically) simpler bound- 
ary conditions, it permits us to use exact boundary conditions rather than the asymptotic 
conditions that are needed when an arbitrarily truncated computational domain is used. 

2. Reduction to a standard eigenvalue problem and the use of a multigrid eigenvalue solver. 

We note that in principle we can solve the quadratic eigenvalue problem as follows |83j . 
Suppose we write equation i|5.9f) schematically as 



(Auj 2 + Buj + C)(t> = 



(5.17) 
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where A, B and C are differential operators. Then we can cast the above equation into a 
standard linear form via the introduction of an auxiliary function, $ defined by 



$ = LUC, 



(5.18) 



Specifically we have 



" A 


" 




' $ " 
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= 0. 



(5.19) 



Hence, our task is now to find an efficient algorithm for solving a standard 2D eigenvalue 
problem [51155] , 

As discussed in Chap.|3 multigrid techniques generally provide the means for very efficient 
solution of multi-dimensional elliptic problems, and with a slight modification of a standard 
multigrid algorithm, we can produce a solver for the eigenvalue problem. A pseudo-code 
of the top-level structure of our multigrid eigenvalue solver for the generalized eigenvalue 
problem 



L(lu, uj 2 )(j> = (Alu 2 + Buj + C) (j) = , 



(5.20) 



is given in Fig. 15.11 Again note that due to the regularity conditions the scalar field must 
vanish on the axis for all k > 0. As a result, for any k, the corresponding family of solutions 
could be naturally parametrized by the modulus of the fc-th order radial derivative of the scalar 
field evaluated at the origin. However, following |22j we have instead chosen to parametrize 
each family of solutions by specifying the value of at some arbitrary fc-dependent point 
(C;s) = (Co(k),0) on the equatorial plane, i.e. (f)(h\ = </>o(Co(^), 0). The other equations of 
System A are solved using a standard multigrid iteration (see Fig. 15.2(1 . We also note that 
we use a FAS (Full Approximation Storage) approach to deal with the nonlinearity of our 
elliptic system. 

A homotopic (continuation) method with initialization by Newtonian solutions: 

As mentioned above, in constructing any particular solution of System A, we will generally 
need a good initial estimate to ensure convergence of the overall iteration (note that Fig. 15. H is 
a pseudo-code for the solution of the equation only, with values of a, f3, ip and a considered 
fixed). We can generate good initial guesses using a simplified version of the homotopy 
(continuation) method [HEj- The idea is to embed the system to be solved (System A) within 
an expanded system. More specifically, if we want to solve a (/-dimensional non-linear system 



F(x) = F(xi, • • • ,x d ) = (Fi(xi,- ■ ■ ,Xd),F 2 (xi, ■ ■ ■ ,x d ), ■ ■ ■ ,F d (xi, ■ ■ ■ ,x d )) = 0, (5.21) 



where the are non-linear functions of (a?i, • • • , x d ) and is the d-dimensional zero vector, 
we consider a family of problems with p parameters, Aj, i = 1 ■ • -p that each take values in a 
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Let 7^ be an initial guess 
repeat 

Solve (4>,L(uj,uj 2 )4>) = for uj 
Rescale so that 4>(( ,0) = (f>r k \ 
Perform Ny multigrid V^-cycles on L(w,w 2 )(/> = 
until \\L(lu,lu 2 )<P\\ < e0 



Figure 5.1: Pseudo-code of the multigrid eigenvalue solver for l|5.9[l . represented here using the 
general form (|5.20() . (4>, L<j>) denotes the inner product of <f> and L<f>. Specifically, 
if the finite difference representations on the mesh (Q,Sj) of <f> and L<fi are (f>i t j and 
(L(f)i : j, respectively, then (4>,L(f>) = J2i j < Pi,j(L(f))ij . <f>^ is the family parameter 
for a particular value of k (see text for additional discussion), is a controllable 
tolerance parameter for the overall algorithm. Ny is another adjustable parameter 
that controls how many multigrid I^-cycles are applied to the discretization of (|5.20() 
at each cycle of the iteration. This iteration is an essential part of the algorithm for 
solving system A, which is itself described in pseudo-code form in Fig. 15.21 

range [aj, b{] for specified constant values a,i and bi. We then replace the above equation with 

G(x;Ai,A 2 ,--- ,A P ) = 0, (5.22) 

where G is defined in obvious analogy with l|5.21[l . such that when A^ = <n for all i, G(x; A,) = 
has a known solution, while when Aj = bi, G(x; Aj) = reduces to the system to be solved, 
namely, F(x) = 0. Therefore, by varying the A, from a, to bi sufficiently slowly 5 and by 
always using the solution from the previous calculation (previous values of Aj) as our initial 
estimate for the new calculation, we can, in principle, obtain the solution to (|5.21(l from the 
known solution to l|5.22(l (i.e. the solution to (|5.22(l for A^ = a^). 

In practice, the continuation parameters, Xi, that we use are N^, N s , 4>(k), m, k and wur, 
where Nq and N s are the number of grid points in the Q and s coordinate directions, re- 
spectively, and wur is the under relaxation parameter defined below. Note that although 
physically k can only take integral values, in the continuation we often find it useful to vary 
it as if it can take on an arbitrary real value. In fact, in some cases we can even start from a 
spherical solution, i.e. with k = 0, and then vary the value of k to get a rotating solution. 

With the continuation method and the other strategies mentioned above, the overall pseudo- 
code for the solution of System A is given in Fig. 15.21 (See the caption for the parameters 
used.) Fig. 15. '6\ shows a schematic representation of the homotopic method in Fig. 15. 21 for two 
typical continuation parameters A^ and Xj . 

5 In general, some of the parameters Ai may not be continuous variable; for instance, one of the parameters might 
be the number of grid points N x in a given coordinate direction, which can only be changed by integral amounts. 
Nevertheless, we have found that variation of such parameters by the admissible small amounts has proven useful 



CHAPTER 5. BOSON STARS IN AXISYMMETRY 



82 



Partition the set {Xi} into subintervals V = {Aj ini | raj = 1 ■ • ■ Ni} 
where An = a, and AjjVi = 

.1=0 

Foreach A = {A lni , A 2 „ 2 , • • • , X, mim } G V 
Call IVP_Solver(A,5 J ) 

3 =j + l 
End 

Subroutine IVP.Solver(A, 5 J ) 
Initialize solutions to 6> J 
I = 
Repeat 

1 = 1 + 1 

Update <fi l and oj according to Fig. 15.11 
Update ip l using standard multigrid 
Update (3 l using standard multigrid 
Update a 1 using standard multigrid 
Update a 1 using standard multigrid 
It 1^0 then 

Call URJJpdate(Z) 
End If 

Until \\u L - u l - 1 \\/\\u l \\ < e , for u € {4>, tp, /3, a, cr} 
Output Solutions to S j+1 
End of Subroutine 

Subroutine UR_Update(7) 
Foreach u 6 {(f), ip, f3, a, a} 

u l = uukU 1 + (1 - ^ur) u 1 ^ 1 
End 

End of Subroutine 



Figure 5.2: Pseudo-code of the initial value solver for axisymmetric rotating boson stars. The 
continuation parameters Xi that are used are N^, N s , , to, k and wxjr- "S* is initial- 
ized cither to a Newtonian solution or to a spherically symmetric solution obtained 
using one of the techniques described in Chap. 0| e is an adjustable convergence pa- 
rameter for the top-level iteration of the algorithm, and is typically set to 10~ 8 . The 
underrelaxation routine URJJpdate, which involves the adjustable underrelaxation 
parameter, wxjr, is used to ensure convergence — see the text for additional details. 
The typical number of overall iterations (loop index I) is between 10 to 50, while the 
number of l^-cycles for each multigrid solver is less than 100. (The actual numbers 
depend on how the partitions of Xi are chosen, i.e. on how good the initial guess S° 
is.) 



in obtaining the final solutions. 
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Figure 5.3: A schematic representation of the nomotopic method shown in pseudo-code form 
in Fig. 15.21 for the case of two continuation parameters, Aj and Aj. The figure 
shows a typical path in the parameter space of A; and Xj. The arrows show the 
stepping-directions through parameter space, while the numbers show the order of 
the stepping sequence. In the first 3 steps we increase the parameter Xi with stepping 
size Aj l7H+1 — Xi. ni (see Fig. I5.2| for the meaning of Aj ini , and note that the stepping 
size is in general different for different i and/or m). At step 4, a further step in the 
Xi direction produces a problem that does not converge, so we then step forward in 
the Xj direction. After step 5 we find that we can again increase Aj until the end of 
step 8, when we realize that increasing either Aj or Xj produces divergence. We thus 
"backtrack" in step 9 and find that we can increase in the Aj direction again. We 
continue with this strategy until we end at step 16, with Aj = &, and Xj = bj. 



Finally, we note that the homotopy method cannot be used without a known initial solution, 
i.e. a solution to 

G(x; oi, a 2 , ■ ■ ■ , a p ) = . (5.23) 

Although we do not have such a solution, we can tune the parameters Aj to the weak field 
limit </>(fc) rs 0, where Newtonian gravity should be a good approximation to general relativity, 
and can then use the solution of the Newtonian problem as the initial guess for 15.23( 1 . 
(Alternately, we can sometimes initialize using a spherically symmetric solution constructed 
using one of the algorithms described in Chap.0J) 

As derived in App. the (non-relativistic) Newtonian limit of the coupled Einstein-Klcin- 
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Gordon system, yields the following coupled Poisson-Schrodinger system: 6 

i<j>(t,x) = -J-V 2 <j) + mV<j), (5.24) 
2m 

V 2 F(t,x) = 47rm 2 00*, (5.25) 
where V is the Newtonian potential. In the stationary case, and assuming an ansatz 

<P(t, x) — > e -' i ( B +™) V( x ) , (5.26) 



the above system reduces to 



Ecf> = -— V 2 + mV0, (5.27) 
2m 



V 2 F = 47rm 2 00*. (5.28) 

Numerical experiments reveal that it is easier to solve ()5.27(l than System A. In fact, using 
an algorithm similar to that described above for System A, we can usually find solutions 
of (|5.27|) using suitably chosen gaussian profiles for </> and V. Once the Newtonian solutions 
(possibly rotating) are in hand, they can be used to start the continuation process for the 
general relativistic solver. 

As a final remark, we mention that the derivation of the Newtonian limit (see App. was 
very useful in identifying which approximations had to be satisfied in the weak field limit. In 
particular, this allowed us to determine whether a particular Newtonian solution was likely 
to be a good initial guess for the full Einstein-Klein-Gordon system. 

4. Underrelaxation: 

Finally, for any value of the continuation parameters, our solution of the 5 coupled elliptic 
initial value equations proceeds via iteration, such that at each cycle (the Repeat loop of 
Fig. I5.2|) each of the 5 unknown grid functions is updated in turn. In the notation of Fig. 15.21 
we have cj) 1 ^ 1 — ► <fr l , tp 1 ^ 1 — > etc. in passing from the I — 1-st iteration to the i-th. 
Empirically, and following Yoshida & Eriguchi [2Hjj we have found it generally useful — and 
in some cases essential — to adopt an underrelaxation strategy, whereby each grid function, u l 
is updated using 

u l = lu vr u 1 + (1 - wur)^ 1 , (5.29) 

where wur is the adjustable underrelaxation parameter and u l is the "bare" solution com- 
puted during the course of the Z-th pass through the main loop of "IVP .solver" . In many 
cases we find that we need wur < 1 for convergence (hence the terminology "underrelax- 
ation" ) , although in some situations wur > 1 can accelerate convergence (as in the usual case 
of the successive overrelaxation technique for the solution of discretized elliptic equations). 



6 As far as we can ascertain, there is no complete derivation of this limit in the literature. In Ref. 1871 Sec. 2.2 an 
attempt was made to derive the Newtonian limit; however, some equations used there are incorrect, and hence the 
derivation is invalid. 
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In general, to achieve optimal convergence, wur must be adjusted on a problem-by-problem 
basis. 

5.2.4 Results 

We now present some results computed using the method described in the previous section for the 
specific cases k = 1 and k = 2. 

Fig. l5.4l shows a typical result for the stationary solution of a boson star with angular momentum 
parameter k = 1. This solution was computed on a finite difference grid with = 129 and 
Nq = 15, and the convergence criterion for the iteration described in Fig. 15.21 was e = 1CP 8 . We 
remind the reader that for k — 1, 2, • • • , the complex field vanishes at the origin, and we thus cannot 
use the central field value to parameterize the continuous family of solutions that exists for each 
value of k. Rather, we parameterize the solutions by choosing (arbitrarily) some point r = ro in 
the equatorial plane (or, equivalently, C = Co)j and then specifying the value of the modulus of the 
scalar field there. Specifically, for the solution shown in Fig. 15. 41 we chose the point £ = 0.5, s = 0, 
and set (j>m = </>o(0.5, 0) = 0.03. The plot at the top of Fig. l5.4l shows the configuration of the scalar 
field, <^o(Ci s )i while the plots below show the various metric components a, (3, tp and a. In all of the 
plots, the axis labeled C corresponds to the equatorial plane, s — 0, or 9 — ir/2, while the far edge 
of each plot, parallel to the £-axis, is s = 1 and is the azimuthal symmetry axis. (As mentioned 
in the introductory section of this chapter, we are restricting attention to equatorially-symmetric 
solutions, and thus only need to solve our equations on a single hemisphere.) The figure shows that 
the matter distribution for the rotating boson star has — in contrast to spherically symmetric boson 
stars or rotating fluid stars — toroidal level surfaces (the toroidal nature of a related configuration 
is clearer in the top plot of Fig. 15.61 which uses the usual cylindrical coordinates (p, z)). We also 
note that the scalar field in this case has a maximum at £ = 0.85, or r = £/(l — C) = 5.7. We 
will see that for k = 2, the maximum value of the scalar field is attained at a larger value of £, 
which implies greater difficulty in constructing solutions at fixed resolution, relative to k = 1. See 
Sec. 15.2. 51 for further details. 

Fig. 15.51 shows a typical stationary boson solution computed for the case k = 2. In this case 
the scalar field and its first radial derivative vanish as r — * 0, and the stellar configuration is 
more compact, in the sense that the support of the modulus of the scalar field on the equatorial 
plane, {£ : |0o(Cj 0)1 > £ } f° r some small number e, is of a smaller interval, with a maximum value 
located further from r = (but still along the equator) than for k = 1. Thus the scalar field at 
(C, s) = (0.5,0) (the point at which we fix <f> for the k = 1 calculations) will be small and, in fact, 
may be comparable to the level of numerical error. In other words, the relative error in </>(2), if fixed 
at (0.5,0), would be large compared to the k — 1 case. Clearly, it is more sensible to choose the 
reference point for 0( 2 ) closer to the point where the maximum of a typical solution in the k — 2 
family is attained. In this case, we have fixed the solution by specifying 0( 2 ) = </>o(0.875, 0) = 0.16. 
(We note that the maximum value, max^ )S |^o(C; s )li °f the configuration shown in Fig. 15.41 is 
attained at C = 0.89, or r = (/(l - C) = 8 - L ) 

From the top plot in the figure, we see that the boson field again takes on a toroidal shape, 
with a more completely evacuated "hole" relative to the k = 1 case. As before, the figures at the 
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Figure 5.4: Typical stationary rotating boson star solution with angular momentum parameter, 
k = 1. The top plot shows the configuration of the scalar field, which has toroidal level 
surfaces. The bottom plots show various metric functions. Note that the solutions 
are symmetric under reflection through the equatorial plane, s — 0. See the text for 
additional details concerning the numerical solutions shown here. In lieu of z-axes 
and labels on the above plots (as well as on the other "surface" plots that appear in 
subsequent figures), we quote the z-extrema for the various data sets visualised here: 

o.o < <MC s) < o.i7, i.o < v(C s ) < i- 1 ) °- 81 < a (C, s) < i-0, o.o < /?(c, s) < 0.014, 

-0.028 < a(C,s) < 0.0. 
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Figure 5.5: Typical stationary rotating boson star solution with angular momentum parameter, 
k = 2. The solutions are symmetric about the equatorial plane. Again, the top 
plot shows the configuration of the scalar field, which has toroidal level surfaces, 
while the bottom plots show various metric functions. As with the k = 1 case, 
the solutions are symmetric under reflection through the equatorial plane, s = 0. 
0.0 < O (C s) < 0.18, 1.0 < V«, s) < 1.2, 0.75 < a((, s) < 1.0, 0.0 < /?(£, s) < 0.014, 
-0.033 < <t(C,s) < 0.0. 
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bottom show the metric functions, which appear qualitatively the same as the k = 1 case except 
that their extrema are located further away from the axis, as we might expect. 

It is also instructive to view the scalar field configurations in ordinary cylindrical coordinates, 
(p,z). Fig. 15. til shows the two scalar field configurations previously displayed, but now plotted in 
(p, z) coordinates. Since the vertical scales of the two plots are very nearly the same, we can see 
from the figure that the scalar field has roughly the same maximum amplitude (<f> max ~ 0.17) in 
the two calculations. It is also clear from the plots that the configuration which ostensibly has the 
larger angular momentum (i.e. the k = 2 calculation, bottom plot) has a larger extent than the 
k = 1 configuration. 

As for the case of spherical symmetric boson stars, we can study some of the properties of our 
rotating configurations as a function of the family parameter (the quantity (f>(u-\ in this case). In 
particular, we can produce plots of total mass vs <j)nA that are completely analogous to Fig. 14.21 
To compute the total mass we use Tolman's expression (modulo our different convention for 
the metric signature): 



MTolman = / (-2T° + T a a ) d 3 X 

ip 6 r 2 sin Qe a 



a 



4> 2 \2uj (uj - (ik) - Tn ot \ dr dO dip 



4n I dr l' 2 d6 ^ r2sin9ea 02 ^ _ /3fc) _ m 2 a 2i 



o Jo 



a 

o 



f 1 f 1 ib 6 e a C 2 
4tt I d( I ds (1 2 c)4 2 [2w (w - (3k) - m 2 a 2 ] . 



Jo Jo 

Note that for any stationary, asymptotically flat spacetime foliated with spacelike hypersurfaces 
that are asymptotically orthogonal to the timelikc Killing field, the Tolman mass expression is 
equivalent to the ADM mass |89l Wi\ . Also note that shows that we have J — kN where J is 
the total angular momentum defined by 

J = J T° 3 ^d 3 x, (5.30) 

and N is the total particle number. Since N has the same qualitative behaviour as a function 
of </>( fe ) as Mroiman (and, in fact the calculations of Yoshida & Eriguchi [2H| — see Fig. 3 of that 
reference — show that the numerical values of N and -Mroiman are quite similar), with maxima of 
both quantities achieved at the same value of 4>(k) > we can at least qualitatively compare the angular 
momenta of stars with varying k using the -Mroiman vs 4>(k) curve and the relation J = kN. 

Fig. 15.71 shows plots of the total mass Mroiman vs the family parameter (j)^ for rotating boson 
stars with k = 1 (top) and k = 2 (bottom). Again, for k = 1 we fix the scalar field value, (f)^ = (p^ 
at (£, s) = (0.5, 0), while for k = 2, where the maximum of (j) tends to be further from the axis, we 
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Figure 5.6: Comparison of scalar field configurations for k — 1 and k — 2. These plots show the 
scalar field solutions previously displayed in Fig. 15. 41 and Fig. 15.51 but now plotted in 
cylindrical coordinates (p, z). The top and bottom plots show the k — 1 and k = 2 
results, respectively. We note that the maximum amplitudes of the two scalar fields 
are comparable for the two cases (<j> m ax ~ 0.17). 0.0 < 4>o(p 7 z) < 0.17 for k = 1, 
0.0 < <f>o(p,z) < 0.18 for k = 2. 
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fix 0(fc) = 0(2) at (£, s) = (0.875, 0). From the plots we can see that in each case, and as for the 
spherically symmetric calculations, there exists a maximum mass limit above which no stationary 
solutions exist. We also note that this mass limit appears to increase with the angular momentum 
parameter, k, which is in accord with naive expectations that centrifugal effects due to rotation 
should partially offset gravitational attraction. 

5.2.5 Some Remarks 

Before ending this section we would like to make a few additional remarks: 

• While Yoshida & Eriguchi's code [2EJ broke down before the maximum-mass star for k = 2 
could be computed, we encounter no difficulties in computing the maximum-mass solution 
(and beyond) in this case (see Fig. I5.7fl . Moreover, we feel that the use of multigrid makes 
our solution strategy potentially more efficient than previous implementations. In particular, 
since the multigrid method is capable of solving elliptic PDEs in O(N) time (where N is the 
number of grid points) , in principle there should be no problems in performing high resolution 
calculations, convergence tests etc. Unfortunately, due to other difficulties we encounter (see 
below), we are currently unable to perform meaningful convergence tests at this time. 

• As discussed in Sec. I5.2.T1 the leading order behaviour for <po(r,0) as r — > is <po(r,8) = 
r k <j) 0k {6) + O{r k+2 ). Therefore MC, %=o = and d^oiC, s)/dCk=o = for i = 1, • • • k- 1. 
In the multigrid routine we only impose the first condition. Similarly for tr(r, 9) the lead- 
ing order behaviour as r — * is <r(r, 0) = r 2 <T2(0) + 0(r 4 ). In other words a(£, s)|c=o = 
and <9cr(C, s)/<9£lc=o = 0. Again we only impose the first condition in the multigrid it- 
eration. Moreover, a closer examination of the r — > limit of System A, performed via 
Taylor series expansion of the variables about r = (for example, we write ip(r,ff) — 
■0(0,0) + r 2 d 2 ip(0, 9)/dr 2 + 0(r 4 ) and similarly for other variables, and then substitute the 
expressions into the governing equations 115.9(1 - 115.1311 1 shows that there are some remaining 
irregular terms (proportional to 1/r) that do not explicitly cancel as r — > 0. Thus it is perhaps 
not surprising that our solutions show irregularity near r = (or cquivalently near £ — 0). 
The fact that we do have some problems with regularity can be seen through close exami- 
nation of the solutions, which reveals non-smoothness near the r = "axis" (see Figs. 15.81 
and I5.9|l . In particular, Fig. 15.91 plots the function da((,s)/d( close to the origin £ = for 
various constant-s slices, s = 0, 0.29, 0.64 and 0.93. The function cr(£, s) was computed in 
a calculation with (pn\ = </>o(0.5, 0) = 0.03, and using a mesh with Nq = 129 and N s = 15. 
Ideally we expect da((, s)/d( — ► as C — > 0. However, the "spikes" in the region < ( < 0.05 
clearly show that the regularity condition is not satisfied. 

• As we have seen in Figs. 15.41 and !5. 51 the larger the value of k, the larger the value of r (or £) 
at which 4>o{Ci s ) is a maximum, and the more concentrated in (£, s) space the stars become. 
Notice that this increased concentration for increased k is largely a function of our use of the 
radially compactified coordinate £. As k increases, so do the overall diameters of the toroidal 
regions, and as r increases the corresponding resolution Ar(£) decreases. Therefore, as k 
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Figure 5.7: The Tolman mass Mxoiman of axisymmetric rotating boson stars with k = 1 and 
k = 2 as a function of family parameter For k = 1 (top) we fix the scalar field 



value. 



(*0 
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at (C,s) = (0.5,0), while for k 



2 (bottom), where the maximum 
2) at (C,s) = (0.875,0). As in 



of <p tends to be further from the axis, we fix <j>i k \ ■ 
the spherically symmetric case, there exists a maximum mass limit above which no 
stationary solutions exist. 
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Figure 5.8: Regularity problem close to the symmetry axis for k = 1. This figure shows a for 
the same solution shown in Figs. 15.41 and 15.61 and plotted in cylindrical coordinates 
(p, z). The figure shows that the function is not completely smooth for small p (i.e. 
close to the symmetry axis on the left). —0.028 < a((, s) < 0.0. 
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Figure 5.9: Irregularity of the metric component cr(£, s) as ( — > for k = 1. Ideally, these curves 
should pass linearly through (0, s) as A£ — * 0. The solution displayed here was 
computed with = <^o(0.5, 0) = 0.03 on a computational domain with = 129 
and N s — 15. The figure shows the derivative of <t(£, s) with respect to £ for constant 
s slices s = 0,0.29,0.64,0.93. The "spikes" in the region < C < 0.05 clearly show 
that the function s) does not satisfy the regularity condition da(0, s)/c?C = 0- 
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increases, higher resolution is needed to sufficiently resolve the stationary configurations. A 
way to remedy this situation is to use a more sophisticated compactification as discussed in 
Sec. 13.51 as well as in the next section. This will increase the complexity of the system of 
equations that need to be solved, but will map the location of the scalar field maximum to 
the central region of the computational domain. In turn, this should partly ameliorate the 
need for increased resolution for larger values of the angular momentum parameter k. 

• The maximum Tolman masses M max = max0 (fc) Mx iman(<^(fc)) we find for rotating boson 
stars are M max = 1.72 for k = 1 and M max = 2.44 for k = 2. In |27| Yoshida & Eriguchi 
found that M max = 1.31 for k = 1 and M max > 2.40 for k = 2. Since their solutions show no 
irregularity at r = 0. the discrepancy (which is rather large) may be due to the regularity 
problems described above. 

• As we have discussed several times, following |27) we parametrize each family of solutions for 
k = 1, 2, • • • by specifying the value of 4>nA at some arbitrary point (£, s) = (Co, 0). Within 
the context of our current particular finite difference solution of System A, this approach has 
the advantage that it is very easy to implement. However, from a theoretical point of view 
it seems more natural to instead specify the A;-th radial derivative of the scalar field. This 
could be implemented by defining a new scalar function (/>o(C, s ) such that [HJ 

MC,s) = <; k MC,s), (5.31) 

and rewriting System A for <po(C) s )- I n this new system we would then specify 9^ <^>o(C; s )lc=o 
as the family parameter, eliminating the ad hoc parameter Co from the algorithm. We discuss 
the initial steps of the implementation of such an approach in the next section. 

• We can replace the individual updates of each variable in the main iteration of our initial data 
solver (the Repeat loop in Fig. I5.2f> by combining the independent multigrid updates shown 
in Figure 15.21 into a single multigrid solve for the five fields. More specifically, in the new 
solver, a basic relaxation step would visit a given grid location (Ci> s j) an d simultaneously 
update the 5 grid function values <fii t3 ■, tpi j ■ ■ ■ Oi j using a 5-dimensional Newton method. 
Thus, each pointwise update in this case would require us to set up and solve a 5 x 5 linear 
system, rather than 5 scalar linear equations. The 5x5 solve is more costly computationally, 
as well as being more tedious to set up since the full Jacobian must be computed. A key 
advantage is that the underrelaxation (UR) part of the algorithm is no longer needed for 
convergence. 

• As in the study of dynamics of boson stars in spherical symmetry, we could study the dynam- 
ics of boson stars with angular momentum in axisymmetry, using the solutions constructed 
in this section as initial data for the dynamical code. However, as mentioned in the introduc- 
tory section of this chapter, the coordinate conditions adopted in graxi (see l|5.4fi|l ) are not 
compatible with our ansatz (|5.7|l . More precisely, if we assume the metric has the form ()5.46|) 
and the boson star has the form l|5.7|l , then the real and imaginary parts of the Klein-Gordon 
equation l|2.45(l will yield different equations for <j)o(r,9). The incompatibility can be traced 
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to the coordinate choice made in adopting the form (|5.46(l whereby the metric component 
9t<p(p, z) is required to vanish |52|. To easily study the dynamics of the rotating boson stars 
constructed here, we would need to adopt a different form for the metric (gt v (p, z) =/= 0), 
which in turn would require substantial modifications to graxi. We have not yet pursued 
this possibility. 

5.2.6 A Proposal for an Improved Rotating Boson Star Solver (System 
B) 

In the previous section we discussed the difficulties we encountered in constructing equilibrium 
solutions of rotating boson stars: we have regularity problems near £ = 0, the maximal values of the 
scalar field are achieved at larger £ on the equatorial plane for increasing k, and the configurations 
are increasingly concentrated in the £-s plane for increasing k. In this section we suggest a way 
of rewriting System A to treat these problems. A new code based on this approach is under 
development, but at this point we do not have any results to report. 
Following j5SJ ED we write the metric in the form 

ds 2 = -e^ +s dt 2 + e 2 ^ (dp 2 + dz 2 ) + e"-y (dip - (idt) 2 , (5.32) 

where 77, S, ip and j3 are functions of the cylindrical coordinates p and z only. (We note that |28U92| 
use spherical-polar coordinates.) The motivation for the above specific form is that a detailed 
analysis (again, by Taylor series expansion of the variables about p = with subsequent substitution 
of the expansions in the equations given below) shows that all equations described below are 
manifestly regular as p — > 0, provided that the p = conditions demanded by local flatness on-axis 
are imposed. Also motivated by regularity considerations, we define a new scalar field variable 
Mp,z) by 

<t>o(p, z) = p k <j> (p,z) , (5.33) 

so that we explicitly factor out the leading order p — > behaviour of the field. The Klein-Gordon 
equation, Hamiltonian constraint, momentum constraint and the equations for 77 and S are derived 
in a manner similar to that used in defining the corresponding equations for System A. For nota- 
tional simplicity, we drop the "0" subscript and the bar so that <j>o(p, z) — > <fi(p, z). We also choose 
k = 1 instead of 87r as previously, and continue to set m = 1. We then have System B: 

^,pp + ^, zz + L P + ^y^j 0,p + il.^.z + (e 2 ^>-\u, + 0k) 2 - e 2 * + + l^^L fA , 

(5.34) 



^,pp + 



— e 

4 



,-2(5 „2 



„2k 



2/c-l 



b 2 p 2k 



(u + Pkf e^- 5 ~ V + ^(1- e 2 ^-'") 



(5.35) 



= 0. 
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/? pp + fizz + (-25 tP + V ,p + ^ 0, P + (V,z ~ 26, z ) >g - 2e^>- r > +s p 2k - 2 <f> 2 k(w + (3k) = , (5.36) 



v.pp + v,zz + v, P 2 + + 7 iz 2 + 4> 2 p 2k 



,2$ _ p 2ip- v -s 



, (5.37) 



S,pp + S tZZ + ( rj.p + - ] 5 p + % z 5, z + 4> p 



2 2fc 



,2ip-ri-8 



-28 



2\ „2 7 ?,P 



(5.38) 



- "(Ap + 2 )p 2 - — = o. 



Further, we compactify both coordinates via coordinate transformations p — > C(p) and z — > //(z) 
given by 



c 



z 

zo + z 



(5.39) 
(5.40) 



The new parameters po and zo can be chosen to map arbitrary coordinate lines p = po and z = zq 
to £ = 1/2 and fi = 1/2, respectively, in the new coordinates (see Sec. 13.50 . Fig. 15.101 illustrates 
the effect of such a transformation. Rewriting System B in (£, /i) coordinates is straightforward 
but leads to lengthy expressions that will not be given here. 

To ensure System B is regular at p = (or equivalently, ( = 0) we apply the locally-flat 
condition 2-0(0, fi) — 7/(0, p) + 5(0, p) — 0. The other regularity conditions are 



0, C (O, p) = 0, C (O, M ) = /3 c (0, /a) - tj )C (0, p) = S x (0, p) = , 
and the other boundary conditions are: 

0(1, M ) = 0(1, /i) = /3(1, p) = r?(l, m) = 5(1, m) = , 



(5.41) 



(5.42) 



<MC, i) = 0(0 1) = £(C, i) = »?(C, i) = <KC, i) = o , 



(5.43) 



<MC, °) = iMC, o) = Ap(C, o) = t/, M (c, o) = ^(c, o) = o . 



(5.44) 



We can now parametrize the family by the value of 0(£, p) = 0o(C A 1 ) at C — and p — 0. The 
multigrid solver will use point-wise simultaneous relaxation of all five unknowns defined at each 
grid point (using a 5-dimensional Newton method) as described in Sec. I5.2."5l 
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Figure 5.10: Demonstration of the effect of adjusting the coordinate-transformation parameters 
p and z . The top figure shows a typical initial configuration of <fio((, m) with p = 1 
and zo = 1. The lower figure shows the same configuration, but now with p = 10 
and z = 10. 0.0 < Mp, z) < 0.01. 
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5.3 Dynamics of Axisymmetric Non-rotating Boson Stars 

In this section we focus on the dynamics of axisymmetric non-rotating boson stars. In the first 
subsection we present the equations of motion that are to be solved, and that have been derived 
using the (2+l)+l formalism outlined in Chap. [2 We note that this derivation is not original to 
this work, nor is the basic code (graxi) that implements the discrete version of the equations 
However, as mentioned previously, we have extended graxi so that it can be initialized with data for 
the complex scalar field representing one or more spherical boson stars in an overall axisymmetric 
configuration. 

In subsequent subsections we perform two classes of numerical experiments. The first class 
involves the head-on collision of equal-mass boson stars, where each star can be boosted towards 
the other via an adjustable initial linear momentum parameter, p z . The second class involves 
perturbation of stable boson stars with a massless real scalar field, generalizing the calculations in 
spherically symmetry that were described in Chap. In both instances, we study the critical be- 
haviour that arises at the threshold of black hole formation. We find evidence for Type I transitions 
in both cases, as well as evidence for scaling of the lifetime, r(p), of near-critical configurations 

r(p) ~ -7 In |p-p*|, (5.45) 

that one normally associates with Type I behaviour (p is, as usual, the family parameter that is 
tuned to generate the critical solution). As far as we are aware, this represents the first instance 
of the observation of Type I critical phenomena in an axisymmetric, general relativistic model. In 
addition, we show evidence for interesting "solitonic" behaviour of interacting boosted boson stars, 
similar to that previously seen in Newtonian computations |32) . 

5.3.1 The equations of motion 

We choose cylindrical coordinates (t,p,z,tp), impose spatial coordinate conditions so that the 2- 
dimensional spatial metric (i.e. the metric in the {p, z) plane) is conformally flat, and adopt 3- 
maximal slicing (i.e. the trace of the 3-dimensional curvature tensor vanishes). As mentioned 
above, we also restrict attention to the non-rotating case. The spacetime metric then has the 
form [U 

ds 2 = {-a 2 + V/ [(/3 P ) 2 + {P Z ?] } dt 2 + 2V> 4 {fi p dp + p z dz) dt + V> 4 (dp 2 + dz 2 + p 2 e 2p °dtp 2 ) , 

(5.46) 

where a, /3 P , (3 Z , ip, and a are functions of p, z and t only. The reason for using a (in contrast 
to the seemingly more natural variable a = pa) as a fundamental metric function is motivated by 
regularity considerations. In particular the leading order behaviour of tr in the limit as p — ► is 
a(p, z) = pai(z) + 0(p 3 ), which can be implemented numerically as a Dirichlet boundary condition. 
On the other hand the leading order behaviour of a is u(p, z) = p 2 o~2{z) + 0(p ), and is much more 
difficult to maintain via finite-differencing of the type used in this thesis. 
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From (|5.46(l we see that the squared-norm, s 2 , of the rotational Killing vector (see l |2.24| l) is 
given by 

s 2 = V/p 2 e 2 " 5 . (5.47) 

Without angular momentum, the twist vector u a (see 12. 25^1 vanishes. 

We write the evolution equation for a in first-order-in-time form by introducing an auxiliary 
variable Ct defined by 



P n = -2^K p p - WK\ (5.48) 

2 „_ f3 z , z - p » 

-n a d a s H '—— 

s 2a 



-n a d a s + — '—z — — — , (5.49) 



where n a is the future-directed, unit-norm, timelike vector orthogonal to the t — const, hypersur- 
faces of the dimensionally reduced 2+1 spacetime, and WK P p , ^K z z are the components of the 
extrinsic curvature on those hypersurfaces. Again, the choice of this particular form for f2 was 
motivated by regularity concerns; certain terms that behave as l/p as p — > 0, and which would 
otherwise appear in the equations of motion for the variable naturally conjugate to a, manifestly 
cancel with the definition 15.48fl . We also note that the leading order behaviour of Ct in the limit 
as p — > is Cl(p, z) = pfti(z) + 0(p 3 ). 

We now adopt a slightly modified form for the stress-energy tensor of the complex scalar field 



TV = T+, = [(V^V^* + V^V^*) - ( V Q 0V Q 0* + m 2 \4>\ 2 + 2A|0| 4 )] . (5.50) 

(Note that the definition differs from that in 12.4411 by an overall factor of 2). In particular, we have 
added a quartic self-interaction term, 2A|</>| 4 to the scalar field Lagrangian, where A is an adjustable 
coupling constant that will satisfy A = for some of the simulations shown below and A = 1 in 
other instances. We note that in addition to the parameter (the particle mass) m which sets the 
length scale of the system, the addition of the self-interaction term introduces a new dimensionless 
parameter, A, into the model, and that our results will generally now depend on the specific value 
of A that is used. In what follows, we have used non- vanishing self-interaction primarily as a vehicle 
to allow us to probe the black hole threshold in certain cases. Additional details will be provided 
in Sec. [5X21 

Wc also define the auxiliary complex scalar variable II via 

Il = n a d a (j), (5.51) 



so that the Klein-Gordon equation for <\> can be recast as a set of first-order-in-time equations. 

In Sec. 15.3.51 we use a real scalar field to perturb the boson stars, analogously to what was 
done in Sec. 14.4.21 in the spherical case. Therefore we again introduce a real massless scalar field 
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03 (t, p, z) with a stress energy tensor 

T$S = 2V^ 3 V^ 3 - . 9ai ,V q 3 V q 3 ■ (5.52) 
The equation of motion for this field is 

V i V^0 3 = O, (5.53) 
and, once more, we define an auxiliary field n 3 

n 3 = n a d a fo , (5.54) 



so that (|5.53|l can be cast in first-order form. 

We note that 0, II, </>3 and II3 are functions of p, z and t only. Since the scalar field (f> = 4>x + *</>2 
and its conjugate momentum II = 11! + in 2 are complex, we now have 12 fundamental real variables 
a, (3 P , (3 z ,^p, ex, Cl, (f) 1 , (f} 2 , 4>3, III, n 2 and II3 that must be evolved. 

As mentioned above (and as for the case of the study of dynamical spherically symmetric boson 
stars in Chap.0J), we choose to implement maximal slicing, so that the trace of the 3-dimensional 
extrinsic curvature tensor of t — const, slices within the 4-dimensional manifold (not to be confused 
with the trace of the 2-dimensional extrinsic curvature tensor ^ K in the 2+1 dimensionally reduced 
spacetime) is zero 

(3) iT = 0. (5.55) 
This condition provides the following elliptic equation for a 

2 (pn, P ) iP 2 + a, zz + a p \2^- + (pa) p j + a z + (pa) 

+8TTil> 4 a (m 2 \(j)\ 2 + A|0| 4 - 2|n| 2 ) - mnaUj = . (5.56) 

The ADM decomposition in the dimensionally reduced 2+1 spacetime will result in one Hamil- 
tonian constraint and two momentum constraints, which, when combined with the demand that the 
(p, z) 2-spaces be conformally flat, provides elliptic equations for i/j, and for (3 P and (3 Z , respectively. 
They are 



+ ^ [W P , P - P z , z ? + (P', P + P p , z ) 2 } + ^ [2apCl + 0P, p - /3%] 2 
= -167T (V 4 (m 2 |0| 2 + AH 4 + |H| 2 ) + |<M 2 + M) - 6 ( P 2 (p*) iP ) p3 

-2 ((p*) tP ) 2 - 2 (pa) tZZ - 2 ({pv), z ) 2 - 16tt (II 2 + 2 , p + , (5.57) 
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j/3'.pp + + -f3\ zp + 32tt—U 3 , p + 16™ (0, p IT + 0* p tf) 

-(fe- 6 T)-<^) < ^ +/ "" ) 



2 /a 

3 \"a 



2ap 



GO 



n, p + 3ft (pa) 



= 0, 



and 



4 1 

/3 2 , PP + -/3^ 2 - ,^P "T 

2ipa) --3{—- 6 — 



2a (p^\ 



Zap f-tj} 



32tt— IT 3 , 2 + 16™ (0, 2 IT + 0* z ll) - 2a(a^ z )p 2 Vt = . 



The evolution equation for a is, with the definition of ft (|5.48(l . given by 

ct = 2/3" {pa) p2 +[3 z a tZ -a(l 



P 



where, as usual, an over-dot denotes a time derivative. The evolution equation for VI is 

1 



Q = 2/3" (pO) p2 + /3^, 2 - — (/3 Z , P 2 - /3", 2 2 ) + 

2a ( A^,P 2 , / - \ 



a ( (V^Xp 



a 



2-0, 2 



+ p<7 2 +(T. ; 



a V 

The definition of II and the Klein-Gordon equation give 

4> = p p (j), p + /3 Z <^ + «n 

and 



+ 647r^p(^) a + 167r^|0 jP | s 





(?) 






a P 









n = /3"n p + /3 2 n, 2 - a<t> (2A|^| 2 + m 2 ) + ^ (?/>>, p + 



i 



a ((/ 9cr ), P <Ap + (P a ), z <!>,*) + §,pp + 4>,zz + — + a, P 0,p + a,z<P,. 



The definition of II3 (|5.54(l and the wave equation for <p 3 l|5.53|) gives 



(5.58) 



(5.59) 



(5.60) 



(5.61) 



(5.62) 



(5.63) 



(5.64) 
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and 



Ha = /5 p n 3 , p + (3 Z U 3 . Z + -n 3 (aptl + 20P p + (5* 



1 



(po-), p </>3, P + (p(r) iZ <fa,. 



The regularity conditions on the symmetry axis are 

a, p (t, 0, z) = V, P (i, 0, z) = {3* p (t, 0, z) = f3 p (t, 0, *) = a(t, 0, z) = fi(t, 0, «) = , 



(5.65) 



(5.66) 



0i, p (t, 0, z) = 02 iP (t, 0, z) = cf> 3jP (t, 0, z) = ni iP (t, 0, z) = n 2lP (t, 0, z) = n 3 , p (i, 0, z) = . (5.67) 
The outer boundary conditions we use, applied at p = p max , z — z max and z = z m j n , are 



a - 1 + 


joa, p + 


za tZ 


= o, 


v> - 1 + 


W\p + 


zip, z 


= o, 


/3 2 + 


pAp + 


z/3% 


= o, 




^ P p + 


zP p z 


= o, 


m ;t + pa 




+ d- 


= o, 


rf2 it + pflp 


+ Z&,z 


+ fl 


= o, 


r^i,t 4- p0i !P - 




+ 4>i 


= o, 


r<f>2,t + p4>2, P - 




+ 02 


= o, 


r<f>3,t + P4>3, P - 


I" Z<t>3,z 


+ 03 


= o, 


rU ht +pU hp + 


zU hz - 


f n x 


= o, 


rH 2 ,t + pU 2 , p + 


^n 2 , 2 - 


f n 2 


= o, 


rn 3 ,t + /9n 3 , p + 


zU 3 ,z - 


f n 3 


= o, 



where r = \J p 1 + z 2 . 

Thus, the system of equations H5.56JI ~ H5.68JI consists of 4 elliptic PDEs (the equations for a, ip, (3 P 
and /3 Z ), and 8 first-order-in time evolution equations (also PDEs) that are hyperbolic in nature. 
The finite differencing of the above equations is given in 0, with trivial modifications for the 
extra terms involving the components 0i,0 2 of the complex scalar fields, and their conjugates 
n 1 ,II 2 . 03 is identified with the original scalar field described in j^j. The finite difference scheme 
also incorporates Kreiss-Oliger style dissipation, with a dissipation coefficient, e<j, similar to that 
defined in IJB.40J) . Again the interested reader can find full details in 8 . 

Following modification to accommodate boson star initial data, the graxi code was subjected to 
convergence tests to ensure that the discrete equations of motion had been implemented correctly. 
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Results from one such test are described in App. [H] 

Since we study dynamics of non-rotating boson stars in this section, we can use static spherically 
symmetric boson star solutions generated using one of the codes described in Chap. 01 as the basis 
for setting up axisymmetric initial data. However, even for the seemingly trivial case where we 
evolve a single spherical star using graxi, the initial data setup is slightly complicated by the fact 
that we need to transform the spherically symmetric solution <fio(r) to the (p, z) coordinates used 
in graxi. We thus generally need to interpolate radially symmetric data onto the discrete domain 
(pi,Zj), possibly translating the data so that the center of the star lies at some specified location 
along the z-axis. Here we use standard polynomial (Lagrange) interpolation which is fourth order in 
the mesh spacings (Ap, Az). Fig. 15. Ill shows a typical result of the interpolation process. Also note 
that all simulations in this section are done with the adaptive mesh refinement (AMR) capability 
provided by graxi, in order to effectively resolve the various length and time scales encountered 
during the dynamical evolution of the scalar fields. In particular, the complex field tends to form 
configurations that are considerable more compact at late times than at the initial time, especially 
in cases where black holes form. 




5.3.2 Head-on Collisions of Boson Stars — Setup of Numerical 
Experiments 

The PDEs solved in the simulations discussed here are those listed in Sec. 15.3.11 The initial data 
for the (first) boson star, always chosen from the stable branch, is given by interpolation of a static 
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spherically symmetric star as just described. The interpolated data, </> Q (p, z), are then scaled by 

$\ p , z )^±$\p,z), (5.69) 

since the definitions 12.44fl and 1)5. 50[) differ by a factor of 2. The star is then "boosted" (roughly 
speaking, given an initial momentum) via (see |87l \'62\ ) 

0( 1 )(p,z;pW)^4 1) (p,z)e^ 1,z , (5.70) 

where p^ 1 is the parameter that controls the magnitude of the boost, and which we will refer to 
as the initial momentum parameter, or, more loosely, as the initial momentum. 

The boosted boson star is then translated so that the center of the star is located at some 
specified point, (0, zi), on the axis of symmetry 

^ 1 \p,z ] p^\z 1 )=T^ 1 \p,z ] p^);z 1 ) , (5.71) 

where T is the translation operator defined by 

T(/( P ,z);z 1 )ee{ /(p ' Z ~ 2l) Z ~ Zl G [Zmin ' Zmax] , (5.72) 
I otherwise 

and z\ is another adjustable parameter. 

Similarly we construct another boosted boson star centered at (0, z%) via 

^ (p, z- pf , z 2 ) ee T (V 2 > (p, z;pW ) ; z 2 ) , (5.73) 

and, finally, the complex scalar field per se is initialized to be the sum of the configurations repre- 
senting the two boosted, translated stars 

0(0, p, z) = 0« (p, z;pW, zi) + (2) (P, z-pf, z 2 ) . (5.74) 

In the simulations described below we have always considered identical stars (i.e. <^>q (p, z) ee 

v o' 



r (p, z) = 4>o{p, z)), each of which is boosted towards the other with the same initial "speed" , (i.e. 
—p^ ee p% ee p z ). We note that the solutions generated from such initial data are symmetric with 
respect to the equatorial plane, but that we do not make explicit use of this fact in the calculations. 
We also note that there is nothing in principle that prevents us from carrying out simulations with 
00 (P> z ) 4>a\p>z) and/or -pi 2) ^ pP . 

In addition to the complex field, <p, we must also supply initial values for the conjugate variable 
II. To this end we define a contribution, IlW(p, z;p^), associated with the first star 

II«(p, z;p«) = ^i0 (1) (p, z;p«) - /3"(p, z)d p ^\p, z;pW) - (3 z (p, z)d z ^\ Pl z;p«) , (5.75) 

where uj\ is the eigenvalue of the first boson star. We similarly define a contribution, Tl^(p, z;pz), 
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for the second star 

n< 2 )(p, z;pW) = Wi^ip, z; P W) - /3"(p, z)d p ^\p, z;p&>) - (3 z (p, z)d z ^ 2 \p, z;p^) , (5.76) 

and then initialize II via 

n(0, /9 ,z)=T(n«(p,z;pi 1 ));z 1 )+r(n( 2 )(p,z;pf );z 2 ) . (5.77) 

The remaining freely specifiable variables cr(t, p, z) and Cl(t, p, z), (which can be associated with 
the gravitational wave content of the spacetime), are initialized to zero: 

<r(0,p,z) = &(0,p,z) = 0. (5.78) 

Finally, initial values for the remaining geometric variables, a,(3 p ,{3 z and ip, are determined from 
the constraint and slicing equations. 

In the simulations we typically use a Courant factor At/Ap = At/Az = 0.2 or 0.3, and the 
dissipation coefficient is = 0.5. 

5.3.3 Head-on Collisions — Solitonic Behaviour 

In the head-on collisions of Newtonian boson stars (i.e. in simulations involving the solution of a 
coupled Poisson-Schrodinger system — see App. |UJ it has been observed that the stars sometimes 
exhibit "solitonic" behaviour (35]. This behaviour is a function of the Newtonian analog of the 
initial momentum parameter, p z , defined above. Specifically, for small initial momenta the stars 
simply merge together to form a single star, while for sufficiently large initial momenta, the stars 
pass through each other as if there was no interaction between them. We have carried out numer- 
ical experiments similar to those described in but using the fully coupled general-relativistic 
equations of motion. We find that solitonic behaviour also occurs in this case for large values of 
the initial momentum. 

In Fig. l5.12l we show results from a typical simulation of a head-on collision of two boosted boson 
stars with no self-interaction (A = 0). The boson stars are initially centered at (0, Zi) = (0,-25) 
and (0, z%) = (0, 25) and have initial momenta —p^ = p^ = 0.4. The computational domain is 
< p < 50, —50 < z < 50. From the figure we see that the stars seem to pass through each other 
as if there were little net interaction between them. 

5.3.4 Head-on Collisions — Critical Behaviour 

We now proceed to study critical behaviour in the context of head-on collisions of boson stars 
with equatorial-plane symmetry. For sufficiently large values of </>o(0) (for instance, (j>o(0) — 0.02), 
variation of the initial momentum parameter, p z , generates a family of solutions that interpolates 
through the black hole threshold; however, in this case the sense of p z is reversed from the normal 
situation — i.e. for small (large) p z black holes do (do not) form. Once we have determined an initial 
bracket for p z that is known to span the threshold, we can, in principle, tune the parameter to the 
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Figure 5.12: Typical solitonic behaviour in a head-on collision of boson stars. The stars are 
initialized via interpolation, translation and boosting of a spherically symmetric 
solution with </>o(0) = 0.02. The stars are initially centered at (0,-25) and (0,25), 
with initial momenta —p^ 1 = p^p = p z = 0.4. The two stars start to overlap at 
t w 31, interfere with each other, and then separate at t ~ 117. Other simulation 
parameters are: At/Az = At/ Ap = 0.2, Ap = Az = 0.78. Note that the temporal 
spacing between successive snapshots is not constant — the time instants displayed 
have been chosen to illustrate the key features of the simulation. 0.0 < \<j>(t, p,z)\ < 
0.03. 
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critical value p* that generates a critical solution. 

However, in our early tuning experiments — which involved no self-interaction for the complex 
field — we found it difficult to determine whether or not a black hole would form in the evolution 
as p z approached the critical value. 7 The addition of the quartic self-interaction term provides 
the stars with more binding, and reduces the distortion in their shapes during the collision. This 
allows us to identify whether a particular calculation is subcritical or supercritical, and we are then 
generally able to determine a value of p* to machine precision. We note that we only incorporate 
the interaction term into the dynamical equations; i.e. we do not compute static boson stars that 
include the quartic term. 8 Thus, with self-interaction turned on, the initial data prescription 
described above will no longer generate a static solution for a single boson star, and in fact, we 
find that in the evolution of data computed with A ^ the central amplitude of the complex field 
tends to grow in time. However, within the time scale for black hole formation in near-critical 
evolutions (t w 300), and for the value A = 1 that we have adopted for the simulations described 
in this section, this growth in amplitude appears to be insignificant. 

Also, we need to remark at this point that when we speak of the black hole threshold in 
the head-on collision calculations described here, we are referring to what we might call "prompt" 
black hole formation, i.e. whether a black hole forms during the initial merger phase of the collision. 
Calculations which are subcritical with respect to this definition are generated by relatively large 
values of p z , and are characterized by the re-appearance of two star-like configurations following 
the initial merger phase (see Fig. 15.151) . Subsequently, however, the stars will re- merge and then 
typically form a black hole (again, see Fig. I5.1cfy . However, this fact that subcritical simulations 
may eventually form black holes does not interfere with our ability to study the threshold defined 
with respect to prompt collapse. 

We now discuss the results from the study of a specific family of head-on collisions in which 
Type I critical phenomena are observed. As in the calculation discussed in the previous sub-section, 
initial data for any member of the family of solutions describe two identical stars (</>o(0) = 0.02) 
centered at (0,-25) and (0,25), and with initial momenta —p? — pz = Pz- The stars are thus 
well separated from each other at t — 0. The computational domain is —50 < z < 50, < p < 50, 
and the base resolution is (N p , N z ) = (129, 257), yielding base mesh spacings Ap = Az = 0.390625. 
As many as 6 additional levels of 2:1 refinement are employed in the calculations, so that the finest 
available resolution is Ap = Az = 0.006104. The critical parameter value for this family of data is 
p\ « 0.21. 

Fig. 15.131 shows snapshots of \<p(t,p, z)\ from a simulation with p z = 0, so that the stars are 
initially at rest. The time development clearly shows that the two boson stars attract each other 
through their mutual gravitational interaction. At f ~ 400 the stars begin to merge, and subse- 
quently they form a configuration which approximates that of a single star. Following the merger, 
the configuration oscillates — emitting some energy in the form of outgoing scalar radiation in the 
process — then collapses and forms a black hole. 

Fig. l5.14l shows the variation of the maximum value of the modulus of the scalar field, |<^| max (t) = 

7 More specifically, we found it difficult to distinguish between black hole formation associated with the super- 
critical evolutions from black hole formation associated with recollapse of subcritical evolutions. 

8 This would certainly be possible to do — see, for example, 1931 . 
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Figure 5.13: Head-on collision of boson stars withp z = 0. Plotted here is the modulus of the com- 
plex field, \(f>(t, p,z)\. The identical stars (4>o — 0.02) are initially centered at (0, —25) 
and (0,25). They move towards one another due to their mutual gravitational at- 
traction, then start to merge at t as 400. The resulting configuration oscillates and 
radiates energy during the merger phase, and then collapses to form a black hole — an 
apparent horizon is first detected at t = 492 (see Fig. 15 . lfc>|l . See the text for addi- 
tional details of the numerical parameters of the simulation. 0.0 < \4>(t, p, z)\ < 0.19. 
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Figure 5.14: Maximum value of scalar field modulus |</>| m ax(i) = max p 

a head-on boson star collision with p z = 0. |</>| m ax(£) oscillates during the in- fall 
phase (presumably due to gravitationally induced distortions of the stars and/or 
the quartic self- interaction) , then starts to grow rapidly during the merger phase, 
t > 430, signalling the collapse of the configuration (analogous to central density 
increasing in fluid collapse) and the imminent formation of a black hole. 



max PjZ |0(i, p, z)\, vs t, again for the case p z — 0. This function initially has a value of 0.014 (note 
that since l|2.44|) and (|5.50(l differ by a factor of 2, the value here is not 0.02). As the evolution 
proceeds it exhibits oscillations (presumably due to gravitationally induced perturbations of each 
of the stars and/or the quartic self-interaction), until t sa 430, when the stars begin to merge. At 
this time, |</>| m ax(i) begins to rapidly increase, and this behaviour signals the imminent formation 
of a black hole (see Fig. I5.16|) . Fig. 15.151 shows the variation of the ADM mass, -Madm vs t, for 
the same calculation (Madm is computed using (2.67) of [8]). The figure indicates that a small 
amount of mass is radiated from the system during the simulation. 

Finally, Fig. l5.16l shows the time development of the apparent horizon detected in the simulation 
at late times — an apparent horizon is first found at t = 492. We note that this figure simply shows 
the location of the apparent horizon in coordinate space for various times and that much of the 
"dynamics" that is visible is probably attributable to coordinate effects. In particular, the proper 
surface areas of the apparent horizons (which we expect ultimately to be closely tied to the black 
hole mass) are not directly correlated to the coordinate areas covered by each of the contours. 
Computation of the former requires evaluation of an integral involving ip(t, p, z) , and ip is a highly 
dynamical quantity, that tends to rapidly increase in the central regions of the domain at late times. 
Thus, although we have not explicitly computed the area of the apparent horizon, we anticipate 
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Figure 5.15: ADM mass vs time for a head-on boson star collision with p z = 0. Although there 
are slight oscillations in this quantity at early times (most likely related to outer 
boundary effects), the mass shows an overall decrease in time, indicating that energy 
is being radiated from the system. 
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Figure 5.16: Development of the apparent horizon for a head-on collision of boson stars with 
p z = 0. Five snapshots of the coordinate location of the apparent horizon are shown 
for t = 491.8,494.0,496.9,498.8 and 503.2 (an apparent horizon is first detected at 
t = 492). As discussed in the text, much of the "dynamics" of the horizon that 
is visible here is likely attributable to coordinate effects. In particular, contrary 
to what the figure suggests, we expect the area of the apparent horizon to be a 
monotonically increasing quantity. 

that it is non-decreasing. 

Fig. 15.171 shows the time development of \<j>(t, p, z)\ for a marginally supercritical head-on colli- 
sion. Here, the momentum parameter, p z ss 0.21 has been tuned to criticality to within a relative 
precision Ap z /p z — O(10~ 15 ). In this calculation, (and due to the non-zero initial boost) the stars 
merge at a significantly earlier time — t ss 90 — than for p z = (compare with Fig. I5.14|l . Following 
the merger, we estimate that the configuration enters the critical state at t *=s 140, and remains 
in that state (with some oscillations) until a black hole forms. An apparent horizon is detected 
at t = 280, and we note that the final black hole will contain most of the mass-energy of the 
configuration at that time. That is, the smallest-mass black hole that can be formed via variation 
of p z has finite mass, as is characteristic of a Type I transition. 

Similarly, Fig. 15.181 shows the time development of \cj>(t,p,z)\ for a marginally subcritical head- 
on collision. Again, the momentum parameter, p z « 0.21 has been tuned to criticality to within a 
relative precision Ap z /p z — O(10~ 15 ). At early times the evolution here is indistinguishable from 
that of the supercritical data. However, in contrast to the supercritical case, at t « 280 two star-like 
configurations begin to reappear and are clearly distinguishable as separate stars at t w 390. Since 
the initial momenta of the stars is moderate, the overall system is gravitationally bound and the 
two stars re-merge and form a black hole. An apparent horizon is detected at t = 570. 

Finally, Fig. l5.19l shows the time of black hole formation (time of first appearance of an apparent 
horizon) ieH plotted as a function of log \p z — p*\ . The plot provides strong evidence for scaling of 
the lifetime, r, of the critical configuration of the form 



t = -7log|p 2 -p* z \ , 



(5.79) 
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Figure 5.17: Marginally supercritical, head-on boson star collision. The time development of 
\4>(t, Pi z) | is plotted for a supercritical calculation with p z tuned to roughly machine 
precision of the critical value p* ~ 0.21. The stars merge at t « 90, significantly 
earlier than in the simulation shown in Fia. l5.l4l The resulting merged configuration 
remains in the critical state for At « 130, then collapses to form a black hole. An 
apparent horizon is first detected at t = 280. 0.0 < \<f>(t,p,z)\ < 0.13. 
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Figure 5.18: Marginally subcritical, head-on boson star collision. The time development of 
\4>(t,p, z)\ is plotted for a subcritical calculation with p z tuned to roughly ma- 
chine precision of the critical value p* w 0.21. At early times the evolution here 
is indistinguishable from that shown in Fig. 15.171 However, in contrast to the su- 
percritical case, at t w 280 two star-like configurations begin to reappear and are 
clearly distinguishable as separate stars at i « 390. Since the initial momenta 
of the stars is moderate, the overall system is gravitationally bound and the two 
stars re-merge and form a black hole. An apparent horizon is detected at t = 570. 
0.0 < \<p(t,p,z)\ < 0.14. 
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Figure 5.19: Scaling law for near critical evolutions of head-on boson star collisions The time 
of black hole formation (time of first appearance of an apparent horizon), £bh, is 
plotted vs log \p z — p* | . The linearity of the data provides evidence for lifetime 
scaling associated with a Type I transition, as previously seen in the spherically 
symmetric calculations described in Chap. 

with a value 7 = 6.6 computed from a least squares fit. Again, observation of such scaling is 
consistent with a Type I transition in the model. 

There remains the question as how the critical solution observed in this study can itself be 
characterized. Although we have not yet studied this matter in detail, we hypothesize that the 
solution can be described in terms of a static boson star on the unstable branch, superimposed 
with (largely) normal-mode oscillations that are not necessarily spherically symmetric. 

5.3.5 Perturbation of Boson Stars by an Aspherical Real Scalar Field 

In this section we present results from a second study of critical phenomena involving boson stars in 
axisymmetry. In this case the calculations generalize those described in Chap. 0] where a massless 
scalar field was used to induce collapse of the boson star. In particular, the massless scalar field 
configuration is not spherically symmetric in the simulations described below. 

For this class of experiment, the initial data for the complex field, </), is simply a single, stable 
boson star centered at the origin, (p, z) — (0,0). As previously, we use a star with 0o(O) = 0.02. 
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The initial data for the real scalar field takes the form of a "generalized gaussian' 



( 



y/(p- po) 2 + e{z- z ) 2 - R t 
A 



) 



2 



o 



(5.80) 



where A3, p$, zq, A, Rq and e are all adjustable parameters. In the calculations described below 
we set po = 0, zq = 0, A = 4, e = 0.4 and Rq — 25. We note that for po — zq = 0, e controls the 
"ellipticity" of the initial pulse, while Ro controls how far the pulse is from the origin at t = 0. 
The overall amplitude factor, A3 is the parameter that is tuned to generate a critical solution, and 
we find A\ w 0.013. We also note that we specify initial data for the conjugate variable II3 so that 
the massless field is almost purely incoming at the initial time. 

The computational domain is — 75 < z < 75, < p < 75, and the base resolution is (N p , N z ) — 
(193,385), yielding base mesh spacings Ap = Az = 0.390625. We again allow up to 6 additional 
levels of 2:1 refinement, so that the finest available resolution is Ap — Az — 0.006104. 

Fig. 15.201 shows the time development of \<p(t,p,z)\ and </>3 (i, p, z) for a marginally supercritical 
perturbed boson star. For clarity the subplots display only a portion of the computational domain, 
namely —50 < z < 50, < p < 50. The top 6 subplots in the figure depict the evolution of the 
boson star, which approximates the critical solution for a period of At ss 70, before forming a 
black hole. An apparent horizon is first detected at t « 220. The bottom 6 subplots show the 
corresponding evolution of the real scalar field, which leaves the computational domain at t f=a 110, 
at which time the boson star is still in a near-critical configuration. Fig. l5.2j displavs corresponding 
results from a marginally subcritical evolution, in which there are indications that the end state of 
evolution is a stable boson star with large amplitude oscillations. 

Thus, we see that the type of non-spherical perturbation used here leads to the same type of 
behaviour seen in the spherically-symmetric calculations, at least at the qualitative level. Among 
other things, this suggests that the unstable static boson stars do not develop additional unstable 
axisymmetric modes when the restriction to spherical symmetry is relaxed (or, if such modes do 
exist, that their growth rates are quite small in comparison to that of the unstable spherical mode). 

We can again plot the time of black hole formation against log \p — p* | to measure the expected 
lifetime scaling for near critical evolutions, as shown in Fig. 15.221 Once more we verify that the 
expected Type I scaling holds. 

Before concluding this chapter, we want to make two additional observations concerning the 
axisymmetric evolutions that we have performed. The first one concerns non-radial oscillations of 
perturbed boson stars, while the second one concerns the linear momentum of such stars. We note 
that we have not studied either matter in much detail, deferring such work to the future. 

When boson stars are perturbed in a non-spherical manner, we observe that, as expected, the 
stars exhibit oscillations which are not strictly radial. Not surprisingly, non-radial oscillations are 
more pronounced when we use a very aspherical distribution of real scalar field as the perturbing 
agent. In Fig. 15.231 we show the evolution of such a highly non-radial configuration of real scalar 
field. We can see from the figure that the initial data consists of several distinct pulses located 
at different positions in the p-z plane. Specifically, the initial configuration of the massless field is 
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Figure 5.20: Perturbation of a boson star (denoted by an elliptical in-going real scalar field 
(denoted <f>). This figure shows the time development of \<j>(t, p, z)\ and 03 (i, p, z) for 
a marginally supercritical simulation, with tuning parameter A3 ss 0.013 (see (|5.80[l ). 
While the simulation domain is < p < 75, — 75 < z < 75, for clarity only the region 
< p < 50, —50 < z < 50 is shown. After the massless scalar field implodes through 
the boson star, the boson star enters its critical state, oscillates about it for a while, 
and then collapses to form a black hole. An apparent horizon is first detected at 
t « 220. 0.0 < \<j>(t,p,z)\ < 0.13 (denoted as $), -0.056 < cj> s (t,p,z) < 0.021 
(denoted as 0). 
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Figure 5.21: This figure shows the time development of \<f>{t, p, z)\ and 4>s(t, p, z) for a marginally 
subcritical simulation (i.e. for A3 of the order of machine epsilon smaller than the 
value used for the simulation shown in Fig. I5.20"jl . In this case there are indications 
that the late time configuration will be a stable boson star with large amplitude 
oscillations, as seen in the spherical calculations of Chap.0] 0.0 < \<f)(t, p, z)\ < 0.13 
(denoted as $), -0.056 < 4> 3 (t,p,z) < 0.021 (denoted as <j>). 



CHAPTER 5. BOSON STARS IN AXISYMMETRY 



117 




Figure 5.22: Lifetime scaling law for near critical evolutions for boson stars perturbed by non- 
spherical real scalar field. The time of black hole formation £bh is plotted as a 
function of log \p — p* \ . The expected (Type I) scaling of the lifetime of the critical 
configuration is again observed. 
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given by 

<f>3(0,P,z) = A 3 exp 
where 



y/(p - po + pi - p2 - Pz) 2 + e(z - z ) 2 - Rq 



Pi 



r>2 



P3 



A exp 
A exp 
A exp 



A 



cosn<9, (5.81) 



A 



A 



(5.82) 
(5.83) 
(5.84) 



and tan (9 = p/z. We have chosen A = 5, A = 8, z c = 35, n = 8, and have set all other parameters 
to the values used in the previously described set of computations. For clarity, we again display 
only a portion of the solution domain, namely —50 < z < 50 and < p < 50. Fig. 15.241 displays 
the corresponding evolution of the boson star, where now the part of the solution domain given by 
—20 < z < 20 and < p < 20 is shown. The figure plots various contour lines of \<f>(t, p, z)\ as a 
function of time, and non-radial oscillations of the star are clearly visible. 

In Figs. l5.2"5l and l5*.26l we show another evolution of a boson star perturbed by an in-going non- 
radial distribution of real scalar field. The boson star is initially located at (0, —7.8125). Again the 
evolution is solved on a domain — 75 < z < 75, < p < 75, while the figures show a portion of the 
domain, —50 < z < 50, < p < 50. The initial configuration for the real scalar field is the same 
as for the previous calculation, except that we choose n — 9 and restrict the support of the real 
scalar field to z > 0. Thus the real scalar field is not symmetric about the equatorial plane in this 
case. Fig. l5.25l shows that at early times, although the modulus of the complex field oscillates, the 
boson star remains at a fixed coordinate location. At t w 80 most of the real scalar field has passed 
through the boson star, and the boson star apparently starts to move in the positive z direction, 
with the amplitude of the complex field still oscillatory. Fig. l5.26l shows the corresponding evolution 
of the real scalar field. 

Fig. ISTH shows the maximum of the modulus of the complex scalar field, max z \4>(t, 0, z)\, 
on-axis as a function of time t, as well as the z-coordinate of the maximum z^, ax (i) defined by 



4>{t, 0, z£ ax (t)) = max \<f>(t,0,z)\ 



(5.85) 



From the graph we can see that the maximum remains at z as — 8 until < w 80, at which point 
it acquires a non-zero velocity which apparently settles down to a speed of about 0.09. However, 
this apparent motion of the boson star seems to be primarily a coordinate effect. Fig. 15.281 shows 
the value of — min z (3 z (t, 0, z) as a function of time t, as well as the apparent velocity of the boson 
star (which exhibits several jumps since the value max \ <p(t, 0, z)\ is not smooth). The two functions 
appear to be of the same order of magnitude, consistent with the conjecture that the motion is 
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Figure 5.23: Dynamical evolution of non-radial distribution of perturbing real scalar field. The 
initial configuration of the real scalar field is given by (|5.81|l . —0.023 < 03 (t, p, z) < 
0.022. 
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Figure 5.24: Dynamical evolution of boson star as a result of non-radial perturbation by a real 
scalar field. The snapshots show both radial and non-radial oscillations of the star. 
0.0 < \<f>(t,p,z)\ < 0.049. 
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Figure 5.25: Evolution of a boson star as a result of perturbation by a non-radial perturbing 
real scalar field, that is not symmetric about the equatorial plane. The boson star 
appears to remain at the same location z w — 8 until t rts 80, when most of the 
real scalar field has passed through the star. The star then appears to move in the 
positive z direction at a speed of 0.09. 0.0 < \<f>(t, p, z)\ < 0.14 (denoted as <f>). 
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Figure 5.26: Evolution of non-radial, non-equatorially-symmetric distribution of perturbing real 
scalar field. Here the real field, whose support at the initial time is confined to z > 0, 
passes through the boson star at t ~ 80, and leaves the computational domain at 
t w 130. -0.047 < (j> 3 (t,p,z) < 0.090 (denoted as (f>). 
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Figure 5.27: The maximum value of the boson star modulus on-axis max z (|</>(£, 0, z)\) (xlO 2 ) 
and the z location of that maximum as a function of time. The boson star starts to 
move at £ « 80, and attains a constant final speed of 0.09. 



mainly due to the particular coordinates that we have chosen. 
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Figure 5.28: Velocity of boson star and the negative of the minimum of the shift component 
P z (t,0,z) as a function of time t. The two functions are of the same order of 
magnitude, indicating that the movement of the boson star depicted in Fig. 15.271 is 
primarily due to the particular coordinates we have chosen. 



CHAPTER 6. CONCLUSIONS AND FUTURE WORK 



124 



CHAPTER 6 

CONCLUSIONS AND FUTURE WORK 

In this thesis we have performed numerical simulations of general relativistic boson stars in both 
spherically and axially symmetric spacetime. The principal new results are as follows. 

In the spherically symmetric case, we studied Type I critical phenomena associated with boson 
stars. In particular, contrary to some previous claims, we found that the end state of subcritical 
evolution is a stable boson star executing large amplitude oscillations, that can largely be under- 
stood as excitations of the fundamental normal mode of the end-state star. For the particular 
example that we examined in detail, the oscillation frequency of the "post-critical" state was es- 
timated to be a 2 /a 2 « 0.0013, in good agreement with the frequency of the fundamental mode, 
<7q = 0.0014, An overall modulation of the envelope of the normal mode oscillations was observed, 
and the origin of this phenomenon remains unclear. 

In the axisymmetric case we developed an efficient algorithm for constructing the equilibrium 
configurations of rotating boson stars. The algorithm was based on an eigenvalue multigrid method 
combined with homotopic/continuation techniques. Complete families of solutions for angular 
momentum parameter values k = 1 and k = 2 were found. The maximum mass computed for the 
two families were l.TM^Jm and 2AM 2 x /m, respectively. The value for k = 1 is only in moderate 
agreement (at the 25% level) with previous work by Yoshida & Eriguchi, who quoted 1.314Mpj/m. 
We are unsure of the origin of the discrepancy at this time, although it may be at least partly 
attributable to the problems with regularity that we encountered with our current code. The value 
for k — 2 is a new result. 

We also demonstrated the existence of Type I critical phenomena in axisymmetry using two 
classes of computations involving non-rotating boson stars. The first involved the head-on colli- 
sion of boson stars, while the second used perturbations of the boson star from a non-spherical 
distribution of massless scalar field. In both cases we measured lifetime scaling from near-critical 
evolutions of the type expected for Type I collapse. 

We also displayed some typical non-radial oscillations of boson stars induced by the gravitational 
interaction with non-spherical distribution of massless scalar waves. 

As for the future, there are several directions in which we can extend our work. For the case 
of rotating boson stars, more work first needs to be done to properly implement the regularity 
conditions. Combined with a modification of the system of equations as suggested in Sec. I3.5l that 
will map the maxima of the unknown functions to the interior of the computational domain, this 
should allow us to find families of solutions with higher values of the angular momentum parameter. 
More importantly, all of our numerical evolutions of boson stars are currently restricted to the non- 
rotating case. It will be very interesting to see how angular momentum affects the behaviour of 
the system. On the other hand, in this thesis we have focused on the dynamical evolution of 
axisymmetric boson stars within the context of critical phenomena. More work can be done to 
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understand the non-radial oscillatory modes of boson stars, and to compare with the corresponding 
modes for neutron stars or white dwarfs. 



BIBLIOGRAPHY 



126 



BIBLIOGRAPHY 



[1] S. H. Hawley. Scalar Analogues of Compact Astrophysical Systems. PhD thesis, The University 
of Texas at Austin, (2000). 

[2] S. H. Hawley and M. W. Choptuik. Boson stars driven to the brink of black hole formation. 
Phys. Rev., D62T04024, (2000). 

[3] D. J. Kaup. Klcin-gordon gcon. Phys. Rev., 172:1331, (1968). 

[4] P. Jctzcr. Boson stars. Phys. Rep., 220:163, (1992). 

[5] L. Lchncr. Numerical relativity: A review. Class. Quantum Grav., 18:R25-R86, (2001). 

[6] M. J. Berger and J. Oligcr. Adaptive mesh refinement for hyperbolic partial differential 
equations. J. Comp. Phys., 53:484, (1984). 

[7] M. W. Choptuik. Universality and scaling in gravitational collapse of a massless scalar field. 
Phys. Rev. Lett., 70:9, (1993). 

[8] F. Prctorius. Numerical Simulations of Gravitational Collapse. PhD thesis, The University of 
British Columbia, (2002). 

[9] M. W. Choptuik, E. W. Hirschmann, S. L. Liebling, and F. Prctorius. An axisymmetric 
gravitational collapse code. Class. Quantum Grav., 20:1857-1878, (2003). 

[10] F. E. Schunck and E. W. Mielke. General relativistic boson stars. Class. Quantum Grav., 
2O:R301-R356, (2003). 

[11] T.W. Baumgarte and S.L. Shapiro. Numerical relativity and compact binaries. Phys. Rep., 
376:41-131, (2003). 

[12] J. A. Wheeler. Geons. Phys. Rev., 97:511, (1955). 

[13] R. Ruffini and S. Bonazzola. Systems of self gravitating particles in general relativity and the 
concept of an equation of state. Phys. Rev., 187:1767, (1969). 

[14] B. K. Harrison, K. S. Thorne M. Wakano, and J. A. Wheeler. Gravitation Theory and Gravi- 
tational Collapse. University of Chicago Press, Chicago, (1965). 

[15] S. Weinberg. Gravitation and Cosmology : Principles and Applications of the General Theory 
of Relativity. Wiley, (1972). 

[16] T. D. Lcc and Y. Pang. Stability of mini-boson stars. Nucl. Phys., B315:477, (1989). 



BIBLIOGRAPHY 



127 



[17] M. Gleiser and R. Watkins. Gravitational stability of scalar matter. Nucl. Phys., B319:733, 
(1989). 

[18] E. Seidel and W. M. Suen. Dynamical evolution of boson stars: Perturbing the ground state. 
Phys. Rev., D42:384-403, (1990). 

[19] J. Balakrishna and W. M. Suen. Dynamical evolution of boson stars. II. excited states and 
self-interacting fields. Phys. Rev., D58:104004, (1998). 

[20] S. L. Shapiro and S. A. Teukolsky. Black Holes, White Dwarfs, and Neutron Stars. Wiley, 
(1983). 

[21] Y. Kobayashi, M. Kasai, and T. Futamase. Does a boson star rotate? Phys. Rev., D50:7721, 
(1994). 

[22] J. B. Hartle. Slowly rotating relativistic stars. I. equations of structure. Astrophys. J., 
150:1005, (1967). 

[23] J. B. Hartle and K. S. Thorne. Slowly rotating relativistic stars. II. models for neutron stars 
and supermassive stars. Astrophys. J., 153:807, (1968). 

[24] V. Silveira and C. M. G. de Sousa. Boson star rotation: A Newtonian approximation. 
Phys. Rev., D52:5724, (1995). 

[25] R. Ferrell and M. Gleiser. Gravitational atoms: Gravitational radiation from excited boson 
stars. Phys. Rev., D40:2524, (1989). 

[26] F. E. Schunck and E. W. Mielke. Rotating boson stars. In F. W. Hehl, R. A. Puntigam, 
and H. Ruder, editors, Relativity and Scientific Computing, pages pp. 138-151, Berlin, (1996). 
Springer. 

[27] S. Yoshida and Y. Eriguchi. Rotating boson stars in general relativity Phys. Rev., D56:762- 
771, (1997). 

[28] S. Yoshida and Y. Eriguchi. New static axisymmetric and nonvacuum solutions in general 
relativity: equilibrium solutions of boson stars. Phys. Rev., D55:1994, (1997). 

[29] C. Gundlach. Critical phenomena in gravitational collapse. Phys. Rep., 376:339-405, (2003). 

[30] S. Noble. A Numerical Study of Relativistic Fluid Collapse. PhD thesis, The University of 
Texas at Austin, (2003). 

[31] B. Rousseau. Axisymmetric boson stars in the conformally flat approximation. Master's thesis, 
The University of British Columbia, (2003). 

[32] D. Choi. Numerical Studies of Nonlinear Schrddinger and Klein-Gordon Systems: Techniques 
and Applications. PhD thesis, The University of Texas at Austin, (1998). 

[33] S. W. Hawking and G. F. R. Ellis. The Large Scale Structure of Space-Time. Cambridge 
University Press, Cambridge, (1973). 



BIBLIOGRAPHY 



128 



[34] C. M. Will. Theory and Experiment in Gravitational Physics. Cambridge University Press, 
revised edition, (1993). 

[35] C. M. Will. The confrontation between general relativity and experiment. Living Reviews in 
Relativity, http://relativity.livingreviews.org/Articles/lrr-2001-4, (2001). 

[36] C. M. Will. The confrontation between general relativity and experiment. Astrophysics and 
Space Science, 283:543-552, (2003). 

[37] R. Arnowitt, S. Deser, and C. W. Misner. In L. Witten, editor, Gravitation: An Introduction 
to Current Research. New York, Wiley, (1962). 

[38] J. W. York, Jr. In L. Smarr, editor, Sources of Gravitational Radiation. Seattle, Cambridge 
University Press, (1979). 

[39] K. Maeda, M. Sasaki, T. Nakamura, and S. Miyama. A new formalism of the Einstein equations 
for relativistic rotating systems. Prog. Theor. Phys., 63:719, (1980). 

[40] R. Geroch. A method for generating solutions of Einstein's equations. Jour. Math. Phys., 
Vol.l2(No. 6):918, (1971). 

[41] S. L. Liebling. Nonlinear Field Dynamics in General Relativity: Black Hole Critical Phenom- 
ena and Topological Defects. PhD thesis, The University of Texas at Austin, (1998). 

[42] J. D. Anderson, Jr. Computational Fluid Dynamics: The Basics with Applications. McGraw- 
Hill, (1995). 

[43] H. Kreiss and J. Oligcr. Methods for the approximate solution of time dependent problems. 
Global Atmospheric Research Programme, Publications Series No. 10, (1973). 

[44] M. W. Choptuik. A study of numerical techniques for the initial value problem of general 
relativity. Master's thesis, The University of British Columbia, (1982). 

[45] W. H. Press et al. Numerical recipes in FORTRAN : the art of scientific computing. Cambridge 
University Press, (1992). 

[46] W. Hackbusch and U. Trottenberg. In W. Hackbusch and U. Trottenberg, editors, Multigrid 
methods : proceedings of the conference held at Kln-Porz, November 23-27, 1981. Berlin: 
Springer- Verlag, (1982). 

[47] A. Brandt. Multi-level adaptive solutions to boundary-value problems. Math. Comput., 31:333, 
(1977). 

[48] C. W. Ucberhubcr. Numerical Computation 2: Methods, Software, and Analysis. Springer- 
Verlag, (1997). 

[49] E. Schnetter, S. H. Hawley, and I. Hawke. Evolutions in 3D numerical relativity using fixed 
mesh refinement. Class. Quantum Grav., 21:1465-1488, (2004). 



BIBLIOGRAPHY 



129 



[50] B. Brugmann, W. Tichy, and N. Janscn. Numerical simulation of orbiting black holes. 
Phys. Rev. Lett, 92:211101, (2004). 

[51] M. W. Choptuik, E. W. Hirschmann, S. L. Licbling, and F. Pretorius. Critical collapse of the 
massless scalar field in axisymmetry. Phys. Rev., D68:044007, (2003). 

[52] M. W. Choptuik, E. W. Hirschmann, S. L. Licbling, and F. Pretorius. Critical collapse of a 
complex scalar field with angular momentum, gr-qc/0405101, (2004). 

[53] J. Thornburg. Coordinates and boundary conditions for the general relativistic initial data 
problem. Class. Quantum Grav., 4:1119, (1987). 

[54] E. Scidel and W. M. Suen. Towards a singularity-proof scheme in numerical relativity. 
Phys. Rev. Lett., 69:1845, (1992). 

[55] M. Alcubierre and B. Brugmann. Simple excision of a black hole in 3+1 numerical relativity. 
Phys. Rev., D63:104006, (2001). 

[56] H. Yo, T.W. Baumgarte, and S.L. Shapiro. Numerical testbed for singularity excision in 
moving black hole spacetimes. Phys. Rev., D64:124011, (2001). 

[57] M. Alcubierre et al. Black hole excision for dynamic black holes. Phys. Rev., D64:061501, 
(2001). 

[58] M.A Scheel et al. Toward stable 3D numerical evolutions of black- hole spacetimes. Phys. Rev., 
D66:124005, (2002). 

[59] H. Yo, T.W. Baumgarte, and S.L. Shapiro. Improved numerical stability of stationary black 
hole evolution calculations. Phys. Rev., D66:084026, (2002). 

[60] G. Calabrese et al. Novel finite-differencing techniques for numerical relativity: application to 
black-hole excision. Class. Quantum Grav., 20:L245-L251, (2003). 

[61] T.W. Baumgarte and S.L. Shapiro. Collapse of a magnetized star to a black hole. Astrophys. J., 
585:930-947, (2003). 

[62] D. Shoemaker et al. Moving black holes via singularity excision. Class. Quantum Grav., 
20:3729, (2003). 

[63] G. Calabrese and D. Ncilscn. Spherical excision for moving black holes and summation by 
parts for axisymmetric systems. Phys. Rev., D69:044020, (2004). 

[64] S. H. Hawley. Private Communication (2003). 

[65] S. M. Carroll. Lecture notes on general relativity gr-qc/9712019, (1997). Eq. (7.5). 

[66] M. W. Choptuik. A Study of Numerical Techniques for Radiative Problems in General Rela- 
tivity. PhD thesis, The University of British Columbia, (1986). 



BIBLIOGRAPHY 



130 



[67] C. W. Misner, K. S. Thorne, and J. A. Wheeler. Gravitation. W.H. Freeman and Company, 
New York, (1973). 

[68] L.I. Pctrich, S.L. Shapiro, and S.A. Teukolsky. Oppenheimer- Snyder collapse with maximal 
time slicing and isotropic coordinates. Phys. Rev., D31:2459, (1985). 

[69] M.W. Choptuik. Consistency of finite-difference solutions of Einstein's equations. Phys. Rev., 
D44:3124, (1991). 

[70] Evans C. R. A method for numerical relativity: Simulation of axisymmetric gravitational 
collapse and gravitational radiation generation. PhD thesis, The University of Texas at Austin, 
(1984). 

[71] R. Friedberg, T.D. Lee, and Y. Pang. Mini-soliton stars. Phys. Rev., D35:3640-57, (1987). 

[72] J. Ventrella. A Numerical Treatment of Spin-\ Fields Coupled to Gravity. PhD thesis, The 
University of Texas at Austin, (2002). 

[73] R. D'Inverno. Introducing Einstein's Relativity. Oxford University Press, New York, (1992). 

[74] A. P. Lightman et al. Problem Book in Relativity and Gravitation. Princeton University Press, 
(1975). Problem 15.3. 

[75] D. Choi, http://lheawww.gsfc.nasa.gov/~choi/axi/axi.html. 

[76] J.L. Friedman, J.R. Ipser, and L. Parker. Models of rapidly rotating neutron stars. Nature, 
312:255, (1984). 

[77] MA. Ruderman and P.G. Sutherland. Rotating supcrfluid in neutron stars. Astrophys. J., 
190:137-9, (1974). 

[78] B.J. Owen et al. Gravitational waves from hot young rapidly rotating neutron stars. Phys. Rev., 
D58:084020, (1998). 

[79] S.I. Yoshida, S. Karino, S. Yoshida, and Y. Eriguchi. A numerical study of the r-mode insta- 
bility of rapidly rotating nascent neutron stars. MNRAS, 316:Ll-4, (2000). 

[80] T.E. Strohmayer. Gravitational waves from rotating neutron stars and evaluation of fast chirp 
transform techniques. Class. Quantum Grav., 19:1321-1326, (2002). 

[81] F.E. Schunck and E.W. Miclkc. Rotating boson star as an effective mass torus in general 
relativity. Phys. Lett. A, 249:389-94, (1998). 

[82] S. Chandrasekhar. The Mathematical Theory of Black Holes. Oxford University Press, (1983). 

[83] L. Komzsik. Implicit computational solution of generalized quadratic eigenvalue problems. 
Finite Elements in Analysis and Design, 37(10):799-810, (2001). 



BIBLIOGRAPHY 



131 



[84] S.R. Kohn and S.B. Baden. The parallelization of an adaptive multigrid eigenvalue solver with 
LPARX. In Proceedings of the Seventh SIAM Conference on Parallel Processing for Scientific 
Computing, page 552, (1995). 

[85] S. Costiner and S. Ta'asan. Adaptive multigrid techniques for large-scale eigenvalue problems: 
Solutions of the Schrodinger problem in two and three dimensions. Phys. Rev., E51:3704-3717, 
(1995). 

[86] R. Burden and J. Faires. Numerical Analysis. Brooks-Cole Publishing, 7th edition, (2001). 
pp.635. 

[87] R. Guenther. A Numerical Study of The Time Dependent Schodinger Equation Coupled with 
Newtonian Gravity. PhD thesis, The University of Texas at Austin, (1995). 

[88] R. C. Tolman. On the use of the energy-momentum principle in general relativity. Phys. Rev., 
35:875, (1930). 

[89] R. M. Wald. General Relativity. University of Chicago Press, Chicago IL, (1984). 

[90] A. Ashtekar and A. Magnon-Ashtekar. On conserved quantities in general relativity. 
J. Math. Phys., 20:793-800, (1979). 

[91] M. Choptuik. Private Communication (2004). 

[92] F. D. Ryan. Spinning boson stars with large self-interaction. Phys. Rev., D55:6081-6091, 
(1997). 

[93] M. Colpi, S. L. Shapiro, and I. Wasserman. Boson stars: Gravitational equilibria of self- 
interacting scalar fields. Phys. Rev. Lett, 57:2485-2488, (1986). 

[94] http: / /www. netlib.org/ eispack. 

[95] B. F. Schutz. A First Course in General Relativity. Cambridge University Press, (1985). 
pp.200. 

[96] F. Pretorius, Unpublished work, (2004). 



APPENDIX A. FINITE DIFFERENCE OPERATORS 



132 



APPENDIX A 

Finite Difference Operators 



Here we list the definitions of all finite difference operators in ID that are used in the thesis (the 
finite difference operators in 2D are defined in a similar way): 
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Table A.l: Finite difference operators used in the thesis. 



We also define ii ± which has the same definition as jj, ± but which has a higher precedence than 
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other algebraic operations, e.g., 

_ r (fix (M r + /;)(M r + 
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APPENDIX B 

Equations of Motion and Finite Difference 
Formulae for the Spherically Symmetric 
Einstein-Klein-Gordon System in Maximal- 
Isotropic Coordinates 



In this appendix we summarize the system of equations governing the spherically symmetric 
Einstcin-Klcin-Gordon model in maximal-isotropic coordinates, as well as the finite difference ap- 
proximations (FDAs) that are used to discretize the system. 
In maximal-isotropic coordinates the spacetime metric is 

ds 2 =(-a 2 + ^p 2 )dt 2 + 2^pdtdr + ^(dr 2 +r 2 dn 2 ) . (B.l) 

We also define auxiliary fields, $ 4 and H 

^ = ti, (B.2) 

n, = ^ (J>i - P4>') , (B.3) 

where i = 1, 2 correspond to the real and imaginary parts of the massive complex scalar field (boson 
stars) , while i = 3 corresponds to the massless real scalar field that is used to gravitationally perturb 
boson stars in our study of critical collapse of those stars. The Hamiltonian constraint, momentum 
constraints, Klein-Gordon equations, and the maximal-isotropic coordinate conditions give 
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rip 2 dr 2 ^ ^ ^ 



The regularity conditions are 



47 rm V £ 4>l - Stt £ II 2 - |(V 2 ^) 2 a = , 



The outer boundary conditions are 
lim ip(t,r) = 

r— >oo 

lim a(t, r) = 

1 — >oo 

lim (3(t, r) = 



i=l i=l 
ri'fY =§alT r . 



^(t,0) = 0, 

iT r (t,0) = 0, 

a'(i,0) = 0, 

#(t,0) = 0, 

n^,o) = o. 
i + £W +0(r -2 ); 

r 



lim — - 



1 = 1 



2C(t) 



+ 0(r- 2 ), 



D(t) 



+ 0(r- 2 ), 



(B.9) 
(B.10) 

(B.11) 
(B.12) 
(B.13) 
(B.14) 
(B.15) 

(B.16) 
(B.17) 
(B.18) 



+ * +— - 0, 
r 
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for some functions C(i) and -D(t). The mass aspect function is defined by 

M(t, r) = {^pj K\ 2 - 2 A 2 (V> + r^') , 
and the location of the apparent horizon is given by that value of r (if any) such that 

4rV>' + 2V> + K r r %i?r = o . 
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In addition, the following evolution equations are used for setting inner boundary conditions when 
black hole excision is used: 
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The finite difference approximations of the Hamiltonian constraint, momentum constraints, 
Klein-Gordon equations, and the maximal-isotropic coordinate choice are given by 
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where r. t = (rj + rj-\)/2. 

The regularity conditions are implemented as 
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for all i and n. The outer boundary conditions are 
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We also adopt a scheme for numerical dissipation given by Kreiss and Oliger 03] ■ I n other words 
an additional term 

(a 4 <0 



^ 4 



is added to the right hand side of (|B.28(1 . for 3 < j < N r — 2 (and similarly for l|B.29|) for Hi), 

where A is defined by 

+ko J 

t n £j / n n n n n \ 

A it = — ( u -4u +6u -4m +m . (B.40) 

Here, e ( i is an adjustable parameter satisfying < < 1, and is typically chosen to be 0.5. 
We note that the addition of Kreiss-Oliger dissipation changes the truncation error of the FDAs 
at 0(At 3 , Ar 3 ) and thus does not effect the leading order error of a second order (0(At 2 , Ar 2 )) 
scheme. The dissipation is useful for damping high frequency solution components that are often 
associated with numerical instability. 
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APPENDIX C 

Convergence Test of the Spherically Sym- 
metric Evolution Code 

Here we present the results of a convergence test of the code that evolves boson stars in spherical 
symmetry. (See Chap. 01 for results computed using this code, and App. [E] f° r full details of the 
equations of motion and finite difference approximations used.) 

The ADM mass is one of the quantities which is most useful for diagnostic purposes. In Fig. IC.ll 
we plot the mass aspect function at the outer boundary of the computational domain, M(t, r max ), 
as a function of time, and from four computations with grid spacings, Ar, in a 8:4:2:1 ratio. The 
simulation involves a pulse of massless scalar field imploding onto a stable boson star as in the 
calculations of the black hole threshold described in Sec. 14.41 

The boson star has a central field value, 4>o — 0.01, while the incoming massless scalar field 
pulse is a gaussian of the form (|4. 109(1 with A3 = 0.001, tq = 40 and a = 3. The outer boundary 
is rmax = 300, and N r = 1025, 2049, 4097, 8193. During the time interval 40 < t < 50, the real 
scalar field is concentrated near the the origin and interacts most strongly with the complex field. 
This results in a localized fluctuation of the ADM mass that is evident in the plots. However, 
M(t, r max ) clearly tends to a constant value as the resolution is increased. In addition, from the 
differences of M(t, r max ) computed at different resolutions (e.g. M Ar (i,r max ) — M 2Ar (t 7 r max ), 
M 2Ar (t, r max ) — M 4Ar (t 7 r max ), etc.), we find strong evidence that the overall difference scheme is 
converging in a second order fashion. 
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Figure C.l: Convergence test of the spherically symmetric code. The estimated ADM mass, 
M(t, r max ), is plotted against time, t, for four calculations using numbers of spatial 
grid points, N r , of 1025,2049,4097 and 8193, so that the corresponding mesh spac- 
ings, Ar, are in a 8:4:2:1 ratio. The initial data parameters for the evolution are: 
0o = 0.01 for the complex field, and A 3 = 0.001, r = 40 and a — 3 for the massless 
field (see <|4.109|) '). The mass decreases with time in general, with a significant fluc- 
tuation at 40 < t < 50, when the real scalar field is close to the origin and strongly 
interacts with the boson star. The change in mass tends to vanish as we go to higher 
and higher resolutions. Combining results from the four calculations we find strong 
evidence that the finite difference scheme is second order accurate as expected. 
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APPENDIX D 

Transformation from Areal to Isotropic Spa- 
tial Coordinates 



In this appendix we summarize the transformation from areal to isotropic spatial (radial) coordi- 
nates. This transformation can be used to generate boson stars initial data in isotropic coordinates 
from data computed in areal coordinates. 

Using areal coordinates, the static, spherically-symmetric metric is 

ds 2 = -a(R) 2 dt 2 + a(R) 2 dR 2 + R 2 dtf , (D.l) 



while in isotropic coordinates it is 

ds 2 = -a(r) 2 dt 2 + V>(r) 4 (dr 2 + rW) . (D.2) 
We thus seek the spatial coordinate transformation 



r = r{R) , (D.3) 

that leaves the t, 9 and <p coordinates unchanged. Comparing coefficients of aTl 2 and the two radial 
elements, we have 

R 2 = ijj 4 r 2 , (D.4) 
a 2 dR 2 = ^dr 2 , (D.5) 



which immediately yields 



iP = */-, (D.6) 



r 



Integrating the last equation, we have 

r = Cexp ( ± 



( ± (D.8) 



where C is some constant. The sign and the integration constant can be fixed by choosing — > r 
when R — > oo: 

lim 1 = 1. (D.9) 

R— >oo it 
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It is instructive to see the implication of the above by considering R sufficiently large that the 
scalar field essentially vanishes. In vacuum the Schwarzschild metric (areal coordinates) takes the 
form 

a=il -ir) ■ 



(D.10) 

where M is the total mass inside the sphere of radius R. Expression 1|D.7(I becomes 

dr dR 



(D.ll) 



r y/R 2 - 2MR 
Integrating, we have 

where C is another integration constant. Applying the boundary condition l|D.9|l we have C = 2, 
i.e., 

M N 2 



R = r[l + -) . (D.13) 



Thus, in a vacuum region, we have 
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Therefore, in vacuum, we can write the metric as 



2M\ . „ 2M N 1 



M 

1 + (D.15) 



ds 2 = - ( 1 - — j dt 2 + M - — ) dR 2 + R 2 dn 2 , (D.16) 



in Schwarzschild coordinates, and 

4 



V 1 + M/2r J \ 2r J v 



(D.17) 



in isotropic coordinates. 

From the above expression for the vacuum metric in isotropic coordinates we can rewrite M in 
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terms of tp and hence find 

1 - M/2r 
a ~ 1 + M/2r 
2 



1. (D.18) 



Therefore in isotropic coordinates, the outer boundary condition for a is 

2 

lim a(r) = — — 1 . 

r— >oo tp 

If we have a numerical solution describing a boson star in areal coordinates, we can transform 
the solution to isotropic coordinates by numerically integrating l|D.8|) from the outer boundary 
inwards. In this case we need to express r in terms of quantities known in the areal coordinate 
system. We may achieve this by noting that from (|D.13|) we have 



2 

On the other hand, in Schwarzschild coordinates we have 



R-M + JW^2MR 



2 V a 

Substituting (fD~20|l in (fElll))) we find 



M = ^(l-- \ . (D.2IM 



(D.21) 

a 

In other words, the coordinate transformation from areal to isotropic coordinates is performed by 
solving 



\R=R„ 



l + R 

a 



J B,=R n 



dr 
dR 



r 
a — . 
R 



(D.22) 



On the other hand, under a coordinate transformation 



x = x(x) . 



(D.23) 



we have 



(D.24) 
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Therefore, under the areal to isotropic coordinate transformation, we have 

-a{r f = [gu], = [<? it ] A - ~a(R) 2 , (D.25) 

where the subscripts I and A denote isotropic and areal respectively, and where quantities with 
(without) a tilde are functions of the isotropic (areal) radial coordinate. The above expression 
implies 

a(r(R)) = a{R) . (D.26) 

In other words, the lapse function transforms as a scalar under the specific transformation of areal 
to isotropic radial coordinates. If we now introduce a uniform grid in R 

R, =(i- 1)AR, i = 1, • • • , N + 1 (D.27) 

where AR = R max /N, then after solving the initial value problem in areal coordinates, we have 

on = a{R l ) = &(r(Ri)) = a{n) = on , (D.28) 

where we have defined on, Ti and on by the above relations. Note that on is simply equal to on, but 
that the latter is defined on a non-uniform grid = (whose values are obtained by solving 

the ODE l|D.22|) ). and that therefore we need to do some interpolation if we want to express the 
function a on a grid that is uniform in r. 

We can make an identical argument for the scalar field, <fi, to obtain the scalar field on the 
isotropic coordinate grid. 

For the conformal factor rp we have 

<AM 4 = [<Mi = (f ) 2 l9RR] A - (^) 2 « 2 = (f ) 2 , (D.29) 
and therefore (we drop the tilde since ip is obviously defined in the isotropic coordinate system), 




(D.30) 



as it has to be. Hence on the non- uniform grid, r i; we have 

*^-> = vQlryf' (D - 31) 

and, as before, we then need to do interpolation to obtain ip on a uniform isotropic grid. 
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APPENDIX E 

Numerical Determination of Spherical Bo- 
son Stars Using a Compactified Isotropic 
Coordinate 



As described in Chap. 0] the system of equations that determines spherically symmetric boson 
stars is an eigenvalue problem. In this appendix we describe a numerical solution of this system 
that uses a compactified, isotropic radial coordinate. 

Let <j>(r), ij)(f) and a(r) denote the modulus of the scalar field, the conformal factor, and the 
lapse function respectively. For notational simplicity let $ = <fi', ^ = ip', and A = a' denote the 
derivatives of the above functions with respect to r. The relevant system of equations can then be 
written as 



A 2* _ 

- + — <I> - ^ TO n 



V 2 "0 + 7T 



+ V> 5 ( — + m 2 ) (f> 2 



^2 2* A A , 4 /2w 2 , , > 
V 2 a + —A - 4mp 4 a —5- - m 2 1 - 
tp \ or 



0. 

0, 

0. 



(E.l) 
(E.2) 
(E.3) 



where V = 3-£z(r 2 -^). The regularity/boundary conditions are: 



$(0) 


= 0, 


(E.4) 


*(0) 


= 0, 


(E.5) 


A(0) 


= 0, 


(E.6) 


lim (f)(r) 

r — >oo 


= 0, 


(E.7) 


lim ip(r) 

r — >oo 


= 1, 


(E.8) 


lim a(r) 


= 1 . 


(E.9) 



r — >oo 



We now perform a compactification of the spatial dimension by introducing a new coordinate, 
C, such that 

C = t^, (E.10) 
1 + r 

where £ g [0, 1]. Under the above transformation we have 
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- = (i-C) 2 -, 

dr V ^ dC 

v 2 = (i-c) 4 v 2 , 

where V 2 = 3^ (C 2 ^)- The previous system of equations can be rewritten as 



(E.ll) 

(E.12) 
(E.13) 



A 2* 

a tp 



V> 4 



(i-CY 



tjj 5 



(1-C) 4 \a 



4 \-2+ m ' * 



2^ , „ -0 4 a /2w 2 , , , 
V C a + — A - in 7T—TU \ —- m ) ( t> 



(1 - C) 4 V a 2 



= 0. 



0. 



(E.14) 
(E.15) 
(E.16) 



with similar regularity/boundary conditions except that the r — > oo limit now becomes (, — > 1. A 
second order finite difference form of the above equations is 




= 0, (E.17) 
= 0, (E.18) 
= 0, (E.19) 



forj = 2,--. ,N r -l. 

The discrete regularity/boundary conditions are 




(E.20) 



i Nr = 0, (E.21) 

with similar equations for ib and a . 

3 j 

Equation (|E. 17fl is linear in <j> and is viewed as an eigenvalue equation with eigenvalue — lo 2 and 
eigenvector 0. More specifically, (|E.17Jl can be written as a matrix equation 



i 



(E.22) 
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where Lij is a tridiagonal matrix. The matrix Lij is in general not symmetric. Transformation 
to a symmetric tridiagonal form can be performed using the FIGI routine in the EISPACK 
package, and then (|E.22|I can be solved with the symmetric tridiagonal eigenvalue solver RATQR 
and TINVIT from the same package. Equation l|E.18|) is non-linear in if>j and is solved by a global 
Newton method, with each Newton step involving a tridiagonal linear solve. Finally, (|E.19|I can 
also be solved for ctj using a tridiagonal solver. The sequence of solving l|E.17|) . <|E.18|) and 1E.19|) 
is repeated until the total change in the variables from iteration to iteration is within a prescribed 
tolerance, typically chosen to be e = 10 -10 . 
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APPENDIX F 

Summary of Equations for Stationary, Ax- 
isymmetric Boson Stars and Associated Fi- 
nite Difference Formulae 



In this appendix we give a summary of the system of equations (System A of Chap. [5j that 
determines stationary rotating, axisymmetric boson stars in quasi-isotropic coordinates, as well as 
our finite difference approximation of the system. 
The spacetime metric is written as 

ds 2 = (-a 2 + %jj A r 2 sin 2 9e 2a j3 2 ) dt 2 + ^ (dr 2 + r 2 d9 2 + r 2 sin 2 9e 2a dtp 2 ) + 2t/>V sin 2 0e 2a (3dtdtp 

(F.l) 

and the equations for the five unknown functions a, (3, ip, a and (p (which depend on r and 9 only) 
are 



cot0 



r 2 sin 2 9e 2a 



ot, r 
a 



0£ 2V^A <Pfi_ 

a ip J r 2 



u- (3k 



2ip >r 



(F.2) 



, , i>,ee 
y,rr + —5- 



1 <*,r 1p,r i 

■ 

r a ip 



-2-Ktp 



ip,e 



a 

air 



i\) fi ip 5 r 2 sin 2 9e 2a 
~~2 777? I P,i 



■ia- 



2r 

. 2 



r 
a 



(1 + ra r ) ^ + - (cot 6 + a g) 

a r a 



<i< 4 



r 2 sin 2 9e 2 ° 



(F.3) 



Prr+^f+[ - +3o> - — + 6^f ) /3 r +( 3cot6» + 3CT e - — + 6^f ) Z£ = -16tt0 



a 



4> 



, fc(w - (3k) 

r 2 sin 2 6»e 2<T 
(F.4) 



a.rr H CK.r H jrdt , 



V> 4 r 2 sin 2 6»e 2<T 



cot^ 



-a, a + 



24>,r 



2f 



2q 



A. 



Airaip <fi 2 [ m 



2{u-(3kf 



(F.5) 



0, 



APPENDIX F. SUMMARY OF EQUATIONS FOR 2D BOSON STARS 



148 



CTrr + — T" + °> + 4 — H + - CT r + CT fl + 4 — H - + 2 COt -5- 

H \ 1/) a r / \ V a J r 

f tL + ^L + l\ ^ r+ ft£ + ^£ +cotd \ &\+ — (a. r + cot 0^L) 
v ip a r J \ if) a J r 2 J ar \ r J 



4 



3V>V sin 2 6»e 2CT , 
+— — I A- 



+4tt 



to 2 + ■ 



3fc 2 



ip 4 r 2 sin 2 #e 2cr 
With the coordinate transformation 



P,e 2 



4a 2 



W -[<!>* + 



,.2 



0. 



(F.6) 



C l + r ' 
s = cos , 



(F.7) 
(F.8) 



the equations become 



i-C 



3C 2 (i-Cf(C 2 0,c), C 3 + ((i- S 2 )0. s ), s - 



k 2 



s (l-s 2 )e 2 



+ 



+ 



m 



0, (F.9) 



V> /i-C 



2a V C 



i-C 



3C 2 (1 _ c) 2 ( C 2^ c ) + (1 _ _ ^ 



-a-cr 



7 + (1 - — + / 



V» 5 C 2 (l-s 2 )e 
8^" 



(i-C) 2 (i + C(i-Ckc)«,c + 



(i-C) 2 Ac 2 + 7^( 1 - s2 )A, 2 

((1 - s 2 )(j, s - s) a s 



i-C 



+27rV> 



C 2 (l - s 2 )e 2CT 



(F.10) 
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(i^) 2 {3C 2 (1 - C) 2 (C 2 /3,c), C 3 + [(1 - -s 2 ) P, s } 



+ (i-Cf 



2 - I ' ~ a ,C , >,C 



, (1-C) 3<r, c --*-+6- 



1-C 



-3,s+(l- S 2 ) ( 3( r s -^ + 6% 



Ac 



+ 167T 



fc(^-/3fc) /1-C\ ,2 n 

(I" I ^ =0 ' 



(F.ll) 



1-C 



3C 2 (l-C) 2 (C 2 «c),3 + (l-5 2 )a„ 



+ (W) 4 (^ + *.c|°,c 



+ 



i-C 

2a 



2 r 



-S+(l- S 2 ) 



2^ [ 2^, 



C 2 (i-CrAc 2 + (i-s 2 )A ; 



Airatp 



4 i2 



2 2 (a; - (3kY 
m - 



+ 



= 0, (F.12) 



(^j {3C 2 (1 - C) 2 (C 2 - C ), C 3 + [(1 - * 2 ) vss sa s ] } 



+a-cr 



i-C 

c 



2 r 



(1-**) (a. s +4% + 2^ 



2s 



+ >-o" 



+ 



i-C 

c 



2% / V>,s , a,. 



( i_^) zi + ri 

1 ^ a 



I 2 (1-C) 2 
a C 



(l-C)a,C-^ a , 



+ 3 ^ (1 rf )e2CT [C 2 (i - C) 2 Ac 2 + (i - 



+47T ■ 



m 2 + 



4a 2 
3fc 2 



V> 4 (1 - s 2 )e 2ff V C 



i-C\ 2 (~ si-' 



a 



(i-C) 4 0,c 2 + 



i-C 



2\ i 2 



(1 - s z )ct>, 



(F.13) 
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The regularity conditions for the system at £ = are given by 

4>(o,s) = o, 

^, c (0,s) = 0, 

Ac(M = o, 

a^(0, s) = 0, 

(7(0, s) = 0, 

and boundary conditions at £ = 1 are 

0(1, s) = 0, 

iKM) = i, 

0(1,8) = o, 

a(l,s) = 1, 

£7(1, S) = 0. 

The boundary conditions at s = (the equatorial plane) are 

<MC,o) = o, 

^,.(C,o) = o, 

0AC,o) = o, 

a,«(C,0) = 0, 

MC,o) = o, 

while at s = 1 (the axis of symmetry), we have 

<KC,i) = o, 

^, s (C,i) = o, 

/3, s (C,i) = o, 

Md) = o. 

For the finite difference formulae, since the expressions are lengthy and messy, we simply note that 
we extend the difference operators defined in App.lAlto operators acting in two spatial dimensions 
in the obvious way. More explicitly, we will have 
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, u. . — u 



A u. 



M. . — U. . 



o *,j 2Ai 
etc.. 



The difference equations are then obtained by performing the following substitutions in l|F.9 
(EH: 
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APPENDIX G 

Weak Field Limit: Derivation of the Poisson- 
Schrodinger system from the Einstein-Klein- 
Gordon System 



In this appendix we will consider the weak field limit (Newtonian limit) of the Einstein-Klein- 
Gordon (EKG) system — thereby deriving the Poisson-Schrodinger (PS) system, which represents 
non-relativistic Newtonian boson stars — with all approximations explicitly stated. The main pur- 
pose of this derivation is to provide the background for checking whether the stationary Newtonian 
solutions obtained by solving the PS system provide good initial estimates for the algorithm used 
to determine the stationary general relativistic solutions (see Sec. I5.2.3|) . The assumptions made 
in deriving the Newtonian limit that need to be checked on a solution-by-solution basis will be 
summarized at the end of the appendix. 

Note that here we write the equations for the EKG system as: 

G a/ 3 = nT a p , (G.l) 

V Q V a 0-m 2 = O, (G.2) 
where k is a constant which is conventionally chosen to be 871". 



G.l Poisson Equation 

The derivation of the Poisson equation for the Newtonian gravitational potential from the Einstein 
field equation is quite standard [HS1 E3|> although the right hand side of the equation (matter 
coupling) is usually not worked out in detail, and the order of truncation is often not clearly 
stated. We suppose that spacetime is nearly flat. The metric can then be written as 

9afj = Vafi + h a p , (G.3) 

where 

\Kp\ = 0(e) <1, (G.4) 

and r) a p is the flat metric in Cartesian coordinates r\ = diag (—1, 1, 1, 1). 1 There are two classes of 
coordinate transformations which preserve the above form: the background Lorentz transformations 

1 This simply means that h is a small number when compared to 77. In the context of numerical computation this 
has a well-defined meaning. Also note that we restrict attention to Cartesian coordinates, since in other coordinate 
systems the components, g a p, of the Minkowski metric can be arbitrarily small. For instance, in spherical coordinates, 
9ee = r 2 — * as r — > 0. 
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and gauge transformations |95|. Gauge transformations can be written in the form 

x & = x a + £ a (x f3 ), (G.5) 

where 

101 = 0(e). (G.6) 



The Jacobian matrix associated with (|G.5|) can be written as 

A s fi = S a p + C,p, (G.7) 

with an inverse given by 



A% = S a -C,p + O(e 2 ), (G. 



as is established by the following: 



Now 



= <5% + 0(e 2 ). (G.9) 



= a " ^ ,a)(^ /3 " r ,/9) V + « " ,«.)(«" ~ C V + 0(e 2 ) 

= Va/3 ~ <f ,a VnP - ,0 Van + h a p + 0(e 2 ) 

= V a + h a 0-£ a ,f3-£p, a + O(e 2 ), (G.10) 



h al 3 = h al 3-Z a ,(3-Zp, a + 0(£ 2 ), (G.ll) 

where in the last step we have adopted the notational convention that we lower and raise tensor 
indices using r\ a 0. That is, h a ^ is, by definition, rj a ^rj^ v h fiu . (Note that using this convention we 
also have g a/3 — r] a ^ — h a/3 + 0{e 2 ). In other words, for all quantities of order 0(e), raising and 
lowering of indices using r\ rather than g will lead to deviations of order 0(e 2 ). ) 
By definition, we have 

— -K 0fiv — 1 /3(U,v _ 1 /9w,/j+ i cri/ 1 /3/i ^ 1 071 1 0v , (y*- 1 -*) 

where 

T A fiv = |r? Ap (V,M + W - *W>) + °( e2 ) ■ (G.13) 
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Note that we have used the fact that we work in Cartesian coordinates (and hence that the deriva- 
tives of i] vanish), and that we assume that derivatives of small numbers are also small 

h aPn = 0(e) (G.14) 

(so that we can drop terms such as h af3 hfi 1 .s)- 
Now we have 

R a » v = -\v ap {hp^ v + hp Plllv - - (/i <->!/) + 0(e 2 ) . (G.15) 

The second terms within parentheses in the above cancel, so that we have 

1 



and 



and 



Rafipv — g {h-av/jp. — hp^ap, ~ h a pjju + hp^. av ) + 0{e ) , (G.16) 



Rp v = \{W vJi p- V,M M - ftM M,^ + hpp» v ) + 0{e 2 ) 

= \{h< 1 v ,pp-Uhp v -h i0v + V/,) + 0(e 2 ), (G.17) 



R = (V, - uh ) + 0{e 2 ) , (G.18) 



where we have defined h= W M and Uh = h^^. 
Constructing the Einstein tensor, we then have 

Rfiv - \rip v R = \ (/i" „,/}„ - Uhp v -h,p v + hp,*, " v - Vpvh 7 6, 75 + V^nti) + 0(e 2 ) . (G.19) 
If we define 

h af3 ee h af3 - \p a!i h , (G.20) 

then 

h = -h, (G.21) 

and 

h a(3 = h a(3 -\p ap h. (G.22) 

Now 

Gf3v = ^ (V v ,(ip - ljh,pv - ^h(3 v + ^np„Uh + h.j3u + hpp,* 1 v - - VP^5, 7S 

+ ^ Vl}v ah - r]p v ahj +0(e 2 ) 

= \(-Vhp v + h» v .^ + h^.» u - mv h l5 ^)+0(e 2 ). (G.23) 
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If we now impose the Lorentz (or harmonic) gauge condition 

h' u, iV = 0, (G.24) 

we have 

G af} = ~ah af3 + 0(e 2 ) , (G.25) 

and hence 

Uh af} = -2nT aP + 0(e 2 ) . (G.26) 
The existence of the Lorentz gauge can be shown as follows. From (IG.llfl it follows that 

h & & = h-2C, a +0(e 2 ). (G.27) 

Also, we have 

= ^-^-^-^(ft-2f , 7 ) + 0(e 2 ) (G.28) 

= h a p - £ a ,0 - + VafiC n + 0{e 2 ). (G.29) 

Hence 

- h& _ = ~ ha p ^ _ a ^ Q + 0(g2) ^ (G.30) 
which leads to the inhomogeneous wave equation 

□r = h a(3 ,/3 ■ (G.31) 

Thus the existence of solutions to (|G.31|I implies the viability of the Lorentz gauge condition for 
the "barred" coordinate system, to order 0(e 2 ). (Note that the solutions allow a homogeneous 
term. Hence the Lorentz gauge is in fact a class of gauges. ) 

Since we expect the Newtonian potential V ~ ft-ooj we assume 

\V\=0(e). (G.32) 

Dimensional analysis shows v 2 ~ GM/r = V \v\ = 0(e 1/2 ). We would expect \T ij \ ~ v\T 0l \ ~ 
v 2 \T°°\, or |T 00 | > |T 0l | > \TV\. This roughly gives \h 00 \ > \h 0l \ > \h ij \. Therefore we further 
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rpOi 

yoo = °( £l/2 )' (G.33) 

rpij 

^oo = (G.34) 

loo - °( £l/2 )' (G.35) 

%5 = (G.36) 
Suppose we can ignore all components of h^ v other than h 00 ; then from (|G.26|) we have 

ah 00 = -2np + 0(e 2 ), (G.37) 
where we have defined p = T 00 . Since we are looking for stationary solutions we have 

Dh aP = V 2 h a(3 + 0(e 2 ) . (G.38) 

which will also be true for a slowly changing gravitational field such that /i 00 ,oo/^- 00 ,ij = 0(e). 
Thus we have 

W 2 h 00 = -2/cp + 0(e 2 ). (G.39) 
Comparing with Newton's law of gravity 

W=~p, (G.40) 



we get 



Moreover, 



h w = -W + 0(e 2 ) . (G.41) 

h = -h = -7 l0Q h 00 + 0(e 2 ) 

= h 00 + O(e 2 ), (G.42) 



h m = h m + ^r, 00 h 

= h 0Q - h 00 + 0(e 2 ) 
= lh 00 + O(e 2 ), 
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and 



or 



That is 



Vooh vv + Vuh 11 = h 

_lft00 + 3ft ii = - h w + o{e 2 ) : (G.43) 

h 11 = h 22 = h 33 = h 00 + 0(e 2 ) . (G.44) 

hoo = h u = h 22 = h 33 = -2V + 0(e 2 ) , (G.45) 

or 

V = - 2 SiwV + C»(e 2 ) . (G.46) 

The metric then becomes 2 

ds 2 = -(1 + 2V)dt 2 + (1 - 2V)(dx 2 + dy 2 + dz 2 ) + 0(e 2 )dx> 1 dx 1 ' . (G.47) 
The harmonic gauge implies 

,„ = h" v , v ~ \rf v \ v = + 0{e 2 ) , (G.48) 

or 

-26» V V V = \^ V {-^V U ) + 0(e 2 ) , (G.49) 

i.e., 

Ko = + O(e 2 ), (G.50) 

which is consistent with our assumption of a (near)-stationary solution. In other words, the har- 
monic gauge condition will be satisfied as long as the solutions are (near)-stationary. Also note 
that (jG.14(l becomes 

V a = 0(e). (G.51) 

Moreover, 

mo~h QO + Vuh U =h = -h m + 0(e 2 ) , (G.52) 

which implies 

h u = 0(e 2 ), (G.53) 

and justifies our assumption 1|G.36(I . (Equation (|G.35|) is satisfied for the diagonal metric l|G.47() . 
which implies = 0(e 2 ). ) 

2 Note that if the original metric is written as 

ds 2 = -a 2 dt 2 + ^{dx 2 + dy 2 + dz 2 ) , 
then a 2 = 1 + 2V or a « 1 + V, and i/j 4 = 1 - 2V or tp a 1 - V/2. 
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To evaluate the right hand side of the Einstein equation we consider 



V - 2 K V ^ V ^* + V ^ V ^*) - <V ( V>V Q 0* + m 2 |0| 2 )] . 

(Note that here V denotes the 4-dimensional covariant derivative operator. ) Since 
have o ~ UJ( t > ~ m </>j where we have assumed 



(G.54) 



e~ lu)t , we 



- = l + 0{e 1 ' 2 ) 



(G.55) 



Now 



i ff 00 [2 ff °V; o ^,o - (9 a/S <f>* a <f>,0 + ™ 2 <^)] + 0(\e<t>\ 2 ) 



,00 



i=l 



0(M 2 ) 



(1 - 2V)V*o6o + (1 - 27)(1 + 2V) ^ + m 2 (l - 2V)(p*(j> 



+ 0(M 2 ) 



(1 - W)P fi 4>, + J2 <M,i + ™ 2 (1 - 2V)(t>*ct> 



-o(ie<^r). 



and 



T 



li 



g 11 [2<7 U <^,1 - (g a ^* a ^^ + m 2 0»] + 0(|e 



a=0,2,3 



+ 0(M 2 ). 



Assumptions (|G.33|) . i|G.34|) can be satisfied if 



^- = 0(eV 2 ) 



(G.56) 



or <j> t i/(f> o = 0(e 1 ^ 2 ), since <pfi ~ m</>. Then 

1 



r 



(10 



~ (w 2 + m 2 ) (/>*(/> + O(e0 2 ) 
™V</> + 0(e</> 2 ), 



(G.57) 
(G.58) 



and the Poisson equation becomes 



VV = -m 2 0> + O(e0 2 ). 



(G.59) 



Also note that T u = 0(e(/) 2 ) and T° l = 0(e xl2 4> 2 ), and hence the 00-component of is dominant. 
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G.2 Schrodinger Equation 

To derive the Schrodinger equation for the complex scalar field, we first note that the Klein Gordon 
equation can be written as 



1 



Now 



-g = v/(l + 2V)(l-2V)3 + 0(e 2 ) 
= l-2V + 0(e 2 ). 



The Klein Gordon Equation becomes: 



1 - 2V 



1 - 2V 
l + 2V ( 



'.o 



3 



Thus 



or 



Now let 



.o »=i 



1 - 2V 



2V 



+ O{e 2 4>)=0. 



1-4V 

1 - 2V <r '" u 1-2V 

-(1 - 2V)cj> fi o + (1 + 2V)V 2 (j> - m 2 = 0{e 2 4>) . 



^,x)^$(t jX )e- mt , 
where $(£, x) is to be a slowly varying function of time. Then 

[-(1 - 2V)(*« - m 2 $ - 2im$, t ) + (1 + 2F)V 2 $ - m 2 $] e"* 1 ™' = (9(e 2 $) , 
-(1 - 2V)($« - 2*m$, t ) + (1 + 2^)V 2 $ - 2m 2 F$ = 0(e 2 $) . 



or 



Let us assume 



V 2 $ 

m 2<j> 



0(6 



m<I> 



m 2<j) 



= o( e 2 ). 



Then keeping track of terms up to 0(e), we have 



(G.60) 



(G.61) 



(G.62) 

(G.63) 
(G.64) 
(G.65) 

(G.66) 
(G.67) 

(G.68) 
(G.69) 
(G.70) 



2im$ t + V 2 $ - 2m 2 V<S> = 0(e 2 ) . 



(G.71) 
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or 

!$ f = -^V 2 $+my$ + 0(£ 2 ). (G.72) 
2m 

We now assume 

*(t,x) = $ ( x ) e - lBt , (G.73) 

so that the full scalar field is 

cj) = $ e - imt = <P oe - tmt ~ zEt = $ e"™* , (G.74) 
where w = m + E. Wc then have 

= --^-V 2 <&Q + mV<S> Q + 0(e 2 <5>). (G.75) 
2m 

Note that from lf(T73|) we have $ t - £$. Therefore l|G"69)) and ifUjOl) will be satisfied if 

E = 0(e) , (G.76) 

or 

- = 1 + 0(e) . (G.77) 
m 

Note that (|G.55|) is then automatically satisfied. Also i|G.68|) will be satisfied if 

- 0(e) . (G.78 

G.3 Summary 

The above derivation is lengthy and hence we present a summary of the development. For our 
purposes stated at the beginning of this appendix, we first solve the PS system 



V 2 V 



E$ Q 



for $ and V. 

Then, ignoring 0(e 2 ) terms, we have 



1 9 

-— V^$ + mV$ . 
2m 



(G.79) 
(G.80) 



-4V 









(G.81) 
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or 





-2V 
















-TV 
















-TV 
















-IV 



(G.82) 



The approximate metric for the EKG system is then 

ds 2 = -(1 + 2V)dt 2 + (1 - 2V)(dx 2 + dy 2 + dz 2 ) , (G.83) 



a = l+V, (G.84) 

= (G.85) 

/3 = 0, (G.86) 

<T = 0. (G.87) 



For the scalar field we have 



$0 , (G.; 



to = 7u + E. (G.89) 

Note that we have assumed (Q, £S>J, (IgTH . (ITT~32l . (iG^Hll - ltG^ . JHSHJ), iG"50t . £T5lJ, 
(|G~53|l . IHSSl, (ESS)- CHOI), (IHTSj) and (lUTSIl . In other words, solutions to the PS equations 
do not necessarily give good initial guesses for the EKG system unless the approximations are 
satisfied. Note that for our purpose of finding stationary solutions, not all of the assumptions 
are independent. Specifically, the assumptions that we need to check on a case- by-case basis are 
(RT521) . (RToTTl . iG"5Hl) . KTm and (RTTgl) : 





= 0(6), 


(G.90) 




= 0(6), 


(G.91) 




= 0(6), 


(G.92) 




= 0(6^), 


(G.93) 


V 2 $ 

TO 2 $ 


= 0(6). 


(G.94) 



Also note that this check is, in general, nontrivial to carry out since in (|G.4|) the flat metric is 
expressed in Cartesian coordinates. Thus, we have to transform our solutions — which are com- 
puted in polar coordinates — to Cartesian coordinates in order to verify the validity of the above 
approximations. 
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APPENDIX H 

graxi Unigrid Static Boson Star Code Test 

In this appendix we present the results of a test of the code, graxi, that we use to evolve the 
massive, complex scalar field and massless real scalar field in axisymmetry (see Chap. 01. 

H.l Background 

graxi is based on a unigrid (i.e. no AMR) code described in [H] that solves the Einstein and massless 
scalar wave equations in axisymmetry. AMR was incorporated into that code by Pretorius [HI E] , 
who later extended the model to include a self-interacting complex scalar field [HSj. In original 
work for this thesis, we extended that code so that it could be initialized with interpolated boson 
star initial data as described in Sec. 15.3.21 this produced the final version of the code that we refer 
to as graxi in the thesis. 

Note that the version of graxi described in has been thoroughly convergence tested, so 

we restrict attention to a test that verifies the boson star initial data interface (our current work) 
as well as the modifications for the complex scalar field per se (Pretorius' work). 

H.2 Test Definition 

The test involves evolution of static boson star initial data. Although this problem may sound 
trivial, it is not. In particular, for the complex scalar field it is anything but a time-independent 
calculation, as a glance at the ansatz 14.63(1 

<f>(t,r)=Mr)e- iut , 

will confirm (recall that we evolve the real and imaginary parts of the complex field separately). 
In addition, the initial data, </>o( r i)j that results from solving a discrete form of the spherically 
symmetric Einstein-Klein-Gordon system with the above ansatz, must always be interpolated to 
the cylindrical grid (pj, Z)~) and this naturally introduces a certain level of error (including departure 
from strict spherical symmetry) in the initial conditions 0(0, pj, z^) and 11(0, pj, Zk). This the finite 
difference solutions <fi h are not expected to be precisely time- independent. 

The parameters for the test are as follows. We initialize the complex scalar field to a single 
boson star with </>o(0) = 0.02, centered at the origin (p,z) = (0,0), on the domain < p < 48, 
—48 < z < 48. We use a Courant factor At/Ap — At/Az — 0.3 and a dissipation coefficient 
e d = 0.5. 
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AMR is not enabled for this test. We note that we have also tested the code with AMR 
enabled, but due to our choice of problem parameters (in particular the truncation error threshold 
that controls the overall regridding algorithm) no refinement is triggered and the results are thus 
essentially identical to those presented here. 

In order to perform convergence tests, we use the above initial data and code parameters and 
then perform calculations with N = 64, 128, 256 and 512 (N p = N + 1, N z = 2N + 1). 

H.3 Test Results 

H.3.1 Time-independence of \(f)(t, z)\ and z) 

Figs. IH.ll and lH72l are 'eye-ball" tests in which we display plots of \4>(t, p, z)\ and ip(t,p,z) that 
should be static, and which do exhibit time-independence to "eye-ball" accuracy. The data shown 
here came from the N — 128 calculation. The length of time spanned by the plots corresponds 
to slightly more than 20 light-crossing times, and about 160 oscillations of each of the component 
fields 4>\ and </>2. 

H.3. 2 Time-independence of max pz \(f)(t, p, z)\ 

Fig. EH plots the maximum value of the modulus of the scalar field, max PjZ \<p(t, p, z)\, vs t from 
the N = 128 calculation. The figure shows that the maximum value exhibits periodic oscillations 
as well as a drift that is due to the accumulation of solution error. We fully expect to see such 
oscillations since numerical errors act as a perturbation to the "exact" stationary solution. 1 

It is an instructive exercise to compare the observed period of oscillation with that computed 
from perturbation analysis of the initial boson star as described in Sec. (|4.4.4H . Assuming harmonic 
time-dependence for the perturbed scalar field and metric component 1(4.120(1 . 1(4.121(1 . and assuming 
that the oscillation is in the fundamental mode, we obtain a theoretical value of a 2 sa 0.00035. 
(Time averaging ^ (t, 0) gives (<£o(*,0)) ~ 0.02 x \/4T = 0.071. ) On the other hand, Fig. HO 
shows a period of oscillation T ss 382, corresponding to a frequency a = 2ir/T = 0.0164. The 
average of the lapse function is (a(t,0)) = 0.993, therefore a 2 /a 2 « 0.00032, in good agreement 
with the perturbation analysis value 0.00035. 

H.3. 3 Time-independence of ADM Mass 

Fig. IH.4I shows the ADM mass Madm(0 vs time t for the N — 128 simulation. The general 
increase in mass starting at t « 450 is indicative of an instability, and implies that the evolution 
will generically break down for a sufficiently long evolution. 

^^Note that the notion of "exact" in this case is relative to the discretization used by graxi and not to the 
continuum. Thus, at any resolution, the "exact" stationary solution, should one exist, would be the solution of the 
difference equations for graxi at that resolution, (including all boundary and regularity conditions) with all discrete 
time derivatives set to zero. 
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t = 792 



( = 936 



1 = 1080 






Figure H.l: Unigrid static boson star code test: \<fi(t,p, z)\. The figure shows N = 128 calculation. 

The initial central value of the scalar field is </>o(0) = 0.02. The period of evolution 
shown corresponds to roughly 20-crossing times and 160 oscillations of 4>\ and fo- 
To "eye-ball" accuracy, \4>(t, p, z)\ is time- independent. 0.0 < \4>(t, p,z)\ < 0.02 x -h=. 
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X = 792 t = 936 t = 1080 




Figure H.2: Unigrid static boson star code test: ip(t, p, z). The figure shows N — 128 calculation. 

The initial central value of the scalar field is 4>o(0) = 0.02. As is the case for \4>(t, p, z)\ 
(Fis. IH.lft . ifj(t,p,z) is time-independent to eye-ball accuracy. 1.004 < tp(t,p,z) < 
1.048. 
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t 

Figure H.3: Unigrid static boson star code test: Oscillation of max p2 \<j>(t,p,z)\ as a measure 
of perturbation. Results shown here are from the N — 128 calculation. The initial 
central value of the scalar field is </>o(0) = 0.02. The observed oscillation frequency 
is roughly a 2 rts 0.00032, and is in good agreement with the fundamental mode 
frequency computed from perturbation analysis, a 2 w 0.00035. 
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Figure H.4: Unigrid static boson star code test: Madm(^) for N = 128 calculation. The initial 
central value of the scalar field is (j>o(0) — 0.02. The general increase in mass starting 
at t w 450 is indicative of an instability, and implies that the evolution will gcncrically 
break down for a sufficiently long evolution. In the inset, we show a detail view of 
the initial fluctuations in Madm(£)- 
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Figure H.5: Unigrid static boson star code test: 4-level convergence test of max p , z \<j>{t, p, z) 
The figure shows that max PiZ \<f>(t, p, z)\ converges to a constant as Ap, Az — » 0. 



H.3.4 Convergence tests of max PtZ (\(j)(t, p, z)\) and Madm(^) 



Figs. IHT51 shows the maximum value of the modulus of the scalar field, max PiZ (| <p(t, p, z)\), vs t 
from the calculations at the four different finite difference resolutions, N = 64, 128, 256 and 512. 
Fig. IH.6I similarly shows the ADM mass as a function of time for the four calculations. The two 
figures provide strong evidence that both quantities are converging. Finally, Fig. IH.7I shows the 
convergence Q-factor for Re(0(£, p, z)) = (f>x(t, p, z), where the Q-factor for a general grid function 
u h is defined by 

Q = \\u™-u»\\ 2 ' (H - 1} 

and which has a theoretical value Q = 4 for a second order scheme. Here || • 1 1 2 denotes the I2 norm, 
and h — Ap = Az. 
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Figure H.6: Unigrid static boson star code test: 4-level convergence test of ADM mass Madm(^)- 
The figure shows that Madm(0 generally converges to a constant as Ap, Az — > 0, 
although there are clearly indications of some problems at later times in this plot 
which are almost certainly a result of imperfect boundary conditions. 




Figure H.7: Unigrid static boson star code test: Convergence rate of Re(4>(t, p, z)) = <pi(t, p, z). 

The theoretical value for a second-order convergent scheme is Q — 4. See text for 
further explanation. 



